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Introduction 


Population dynamics [1] constitutes a widespread branch of investigations which 
finds particularly important applications within the realm of life science. In general 
terms, it aims at characterizing the time evolution of interacting species of homologous 
entities, so to unravel, among others, the fundamental mechanisms which drive their 
dynamics as observed in real systems. 

Population is indeed a technical term which is referred to various, completely dis- 
tinct fields of applications ranging, from e.g. the level of expression of a protein in a cell, 
to the number of animals in a finite ecosystem. 

The classical approach to population dynamics [2] relies on characterizing quanti- 
tatively the densities x = (x1,..., n) of n species through a system of ordinary differ- 
ential equation of the type 


d 
ae ASA) 


where the function f depends on the specific interactions being at play. In other words, 
the analytical expression for f, incorporates pure competition, predator-prey interac- 
tions, or even cooperative effects. Moreover, a specific delay might be required to ac- 
count for the processing time which is needed for a system under scrutiny to react to an 
external stimulus or signal [3]. This is a paradigmatic problem of many biological path- 
ways. A classical example is the haematopoiesis, the formation of blood cellular compo- 
nents [4], where the production of thrombocytes is delayed of seven days with respect 
to the level of megakaryocytes, the progenitor cells. These latter take in fact seven days 
to complete their differentiation cycle. Towards a refined level of approximation, more 
than one independent variable is often to be assumed, which in turn implies dealing 
with the partial differential equation for an exhaustive modelization effort [5]. For in- 
stance, when tracing the dispersion of a diffusing chemical compound, space and time 
are to be explicitly encapsulated into the mathematical description. 

However, and despite the degree of coarse-graining intrinsic to the model, all these 
are phenomena that can be tackled via the population viewpoint, namely focusing on 
the evolution of homogeneous family of constituents as whole, and solely allowing for 
effective (global) interactions between microscopic elements. It is customary to refer to 
this level of description as to the deterministic theory. Noise and other disturbances can 
be eventually hypothesized to alter the ideal deterministic, hence reproducible, dynam- 
ics but always act as a macroscopic bias. 

As opposed to this formulation, a different level of modeling can be invoked which 
instead focuses on the individual-based description. This amounts to characterizing 
the microscopic dynamics via explicit rules governing the interactions among individ- 
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uals and with the surrounding environment. This former approach has been recently 
adopted in various contexts such as predator-prey interactions [6], metabolic reac- 
tions [7], and epidemic models [8]. These models are usually implemented numerically 
through algorithms which use random numbers, and for these reason we refer to these 
models with the term stochastic.Th e stochasticity is now intrinsic to the systems and 
stems from the microscopic finiteness of the investigated medium. 

Those alternative, conceptual strategies translate into different tools for character- 
izing a given phenomenon under inspection, and it is therefore of interest to highlight 
similarities, and/or discrepancies, in the associated predictions. A viable method that 
enables us to bridge the gap between the two levels of description is the so-called van 
Kampen’s system size expansion. Starting form the stochastic scenario and perform- 
ing a perturbative development with respect to a small parameter which encodes the 
amplitude of finite size fluctuations, one recovers, at the leading order, the mean-field 
equations. These latter govern the coupled evolution of the average population amount, 
as in the spirit of the deterministic representation. Including the next-to-leading order 
corrections, one obtains a description of the fluctuations, as a set of linear stochastic 
differential equations. Such system can be hence analyzed exactly, so allowing us to 
quantify the differences between the stochastic formulation and its deterministic ana- 
logue. 

Again, let us emphasize that fluctuations do not arise from an external noise. Despite 
the evidence that it is always present in actual population dynamics and that it is an 
essential ingredient of life processes, noise is often omitted. When instead considered, 
it is frequently assumed to act as a source of disorder and it is included in the dynamics as 
an external elements. At variance, the individual-level approach allows us to investigate 
the unavoidable intrinsic noise, which originates from the discreteness of the system and 
that has to be considered in any sensible model of natural phenomena. 

In this thesis we will apply the van Kampen theory to a selection of problems in 
biomedicine and biology, for which we shall also propose dedicated interpretative frame- 
work. More specifically, we will focus on the role of finite-size corrections, and show 
how these might significantly alter the dynamics. In particular we shall analyze the 
molecular mechanisms involved in the perception of pain, and the autocatalytic reac- 
tions which widely occur in many biochemical processes. 

The first chapter is devoted to introducing the main aspects involved in the neurobi- 
ology of pain. We will briefly describe the various types of pain and the corresponding 
pharmacological treatments. Moreover we shall mention the new frontier of the per- 
sonalized medicine, making always reference to the problem of pain emergence. 

The second chapter presents a gallery of methods, which we did developed as sup- 
port of the experimental data analysis. First, a deterministic model is put forward to 
interpret drug kinetics data, with the aim of quantifying the individual response to a 
specific medical treatment. Moreover, the problem of missing data in microarray ex- 
periments is considered, and a method suggested that exploit the similarity between 
sequences. In the last section of this chapter, we report on a new algorithm for detect- 
ing different level of clustering. 

The third chapter deals with a stochastic model to investigate the microscopic pro- 
cesses which trigger the sensation of pain. The model accounts for the action of anal- 
gesic drug and introduces an effect of competition among chemical species populating 
the bloodstream. Regular oscillations in the amount of bound receptors are detected, 
following a resonant amplification of the stochastic component intrinsic to the system. 


Introduction IX 


The condition for such oscillations to occur are studied, resorting to combined numer- 
ical and analytical techniques. Extended and connected patches of the admissible pa- 
rameters space are detected which do correspond to the oscillatory behaviors. These 
findings are discussed with reference to the existing literature on patients response to 
the analgesic treatment. 

In the fourth chapter we present a minimalist stochastic model for the mechanism 
of action of tramadol. The model accounts for the process of metabolization through 
the cytochrome CYP2D6 and the interactions between molecules and target receptors. 
From the master equation, through the van Kampen method, we obtain the macro- 
scopic equations, and the Fokker Planck equation governing the fluctuations around 
the deterministic behavior. The role of fluctuation is discussed with reference to clinical 
tests and outcome of pharmacological therapy. 

The last Chapter is dedicated to a different topic, where similar modelization tech- 
niques prove necessary. A simplified scheme of k coupled autocatalytic reactions is ana- 
lyzed. Th is problem is supposedly related to the inner dynamics of a simplified model of 
cell, as previously recognized in the literature. The role of stochastic fluctuations is elu- 
cidated through the use of the van Kampen system-size expansion and the results com- 
pared with direct stochastic simulations. Regular temporal oscillations are predicted 
to occur for the concentration of the various chemical constituents, with an enhanced 
amplitude resulting from a resonance which is induced by the graininess of the system. 
Space is then accounted for, resulting in organized spatio-temporal structures. 


Chapter 1 
Teory of Pain 


The International Association for the Study of Pain (IASP) defines pain as “an un- 
pleasant sensory and emotional experience which we primarily associate with tissue 
damage or describe in terms of such damage, or both” Although this is a vague def- 
inition, it recognizes that pain is a combined sensory, emotional, and cognitive phe- 
nomenon. Interestingly, physical pathologies do not necessarily manifest when a pa- 
tient experiences pain. Indeed, pain can be thought of as being composed of three hi- 
erarchical levels: a sensory-discriminative component (e.g., location, intensity, qual- 
ity), a motivational-affective component (e.g., depression, anxiety), and a cognitive- 
evaluative component (e.g., thoughts concerning the cause and significance of the pain) 
[9]. Clinically, this conceptual vision is useful to focus the attention on the broad range 
of factors that may contribute to the emergence of pain as reported by the patients. 


In general terms, according to the medical literature, pain is “whatever the patient 
says it is” One may argue that pain is a warning sensation delivered to the patient’s brain 
signalling that a stimulus is causing, or may cause, damage. Following this philosophy, 
the best clinical approach in most circumstances is to assume that the patient is report- 
ing a true experience, even in the absence of an obvious demonstrable source of tissue 
injury. However, accepting a patient’s complaint of pain as valid does not necessarily 
demand the initiation of a specific treatment. 


The difficulties in the pharmacological management of pain are not solely related to 
the complex underlying network of molecular interactions, which remains at present 
to be fully elucidated. As an additional complication, the same administered dose of 
drug can result in a wide variability of effects, depending on specific individual traits. 
This observation points to the need for the so called personalized medicine, one of the 
new challenges of the post genomic era.Th is, relatively novel discipline elaborates in- 
formation on the genetically determined variability in the metabolism of drugs, aiming 
at developing personalized pharmacological protocols to enhance the effectiveness of 
the therapy. 


This chapter is devoted to providing a general, though synthetic, overview on the 
neurobiology of pain. We shall be in particular interested in highlighting the various 
types of pain and discussing possible treatment strategies. The pharmacological aspects 
are here briefly reviewed with reference to the issue of personalized medicine. 
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11 The pathophysiology of pain 


Pain perception, or nociception (from the Latin word for “hurt”), is the process by 
which a painful stimulus is transmitted from the site of stimulation to the central ner- 
vous system. Several successive steps are to be accounted for in the nociception process: 
The contact with stimulus, the reception, the transmission and the pain center percep- 
tion [10, 11, 12]. In the following we shall set down to analyzing those processes. 

The contact with the stimulus can be mechanical (pressure, punctures and cuts) 
or chemical (burns). As an example, we can consider the cut of a finger. The stimu- 
lus is captured by nociceptors, nerve cell endings that, like normal sensory neurons, 
travel in peripheral sensory nerves. Their cell bodies lie in the dorsal root ganglia of 
peripheral nerves just inside the spine. Nociception uses different neural pathways than 
normal perception (like light touch, pressure and temperature). In response to non- 
painful stimulation, the first group of neurons to fire are normal somatic receptors. At 
variance, when external agents cause pain, nociceptors are activated atfi rst. Nocicep- 
tors, indeed, sense pain via free nerve endings rather than specialized endings such as 
those involved in resolving the touch or pressure feelings. However, while normal sen- 
sory neurons are myelinated (insulated) and conduct quickly, nociceptor neurons are 
lightly or non-myelinated and slower. We can divide nociceptors into three classes: 
The Ad-mechanosensitive receptors and the Ad-mechanothermal receptors which are 
faster conducting neurons that respond to mechanical stimuli (pressure, touch) and to 
heat, and the polymodal nociceptors (C fibers) which are slowly conducting neurons 
that respond to a variety of stimuli. After the cut of thefi nger, several factors contribute 
to the reception of pain: The mechanical stimulation, the potassium released from the 
insides of the damaged cells, the prostaglandins, histamines and bradykinin from im- 
mune cells that invade the area during inflammation, and the substance P' from nearby 
nerve fibers. These substances cause action potentials in the nociceptor neurons. The 
first sensation experienced at the moment of the injury is an intense pain. The signal 
for this pain is conducted rapidly by the Ad-type nociceptors. The pain is followed by 
a slower, prolonged, dull ache, which is instead conducted by the slower C-fibers. 

In this way, the signals originated from the injury event travel into the spinal cord 
through the dorsal roots. There, they make synapses” on neurons within the dorsal horn 
(the top half of the butterfly-shaped gray matter). They synapse on neurons within the 
spinal cord segment and also on neurons one to two segments above and below their 
segment of entry. These multiple connections relate to a broad area of the body. This 
in turn explains why it is sometimes difficult to determine the exact location of pain, 
especially internal pain. The secondary neurons send their signals upward through an 
area of the spinal cord’s white matter, termed the spinothalamic tract. This area can be 
pictorially seen as a highway: Traffic from all of the lower segments rides up the spinal 
cord. The signals of the spinothalamic tract is transmitted along the spinal cord through 
the medulla (brain stem) and synapse on neurons which are located in the thalamus, 
and act as brains relay center. Some neurons also synapse in the medulla’s reticular 
formation, a region which is deputed to physical behaviors. Nerves from the thalamus 
then pass the signal to various areas of brain’s somatosensory cortex. 


1 Substance P is a neuropeptide, namely a short-chain polypeptide that functions as a neurotransmitter and 
as a neuromodulator which alters the excitability of the dorsal horn ganglion (pain responsive neurons). 


2 Synapse is the functional connections between neurons, or between neurons and other types of cells. 
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Once the pain information reaches the brain, it gets processed following schemes 
which are still to be completely resolved. Obviously, the signals should eventually reach 
the motor cortex, then down through the spinal cord and to the motor nerves. These 
impulses would cause muscle contractions to move the hand so to escape the source of 
potential injury. 

Physicians and neuroscientists generally classify pain depending on its duration, or 
according to the associated clinical characteristics. 

When it comes to quantifying the time duration, pain can be divided in three main 
classes, namely acute, chronic and cancer pain. Acute pain is caused by an injury to the 
body. It warns on a potential damage that requires immediate action. It can last for a few 
minutes to six months and fade off when the injury heals. Chronic pain persists long 
after the trauma has healed (notice that, in some cases, it may occur in the absence of 
any trauma). Chronic pain does not warn the body to respond, and it usually lasts longer 
than six months. Finally, cancer (or malignant) pain is clearly associated with malignant 
tumors. Tumors tend to invade healthy tissues so exerting a mechanical pressure on the 
nerves or blood vessels, which generates pain. Occasionally, cancer pain can be also 
classified as chronic pain. 

A grouping based on inferred pathophysiology, broadly classifies pain syndromes 
into nociceptive, neuropathic, psychogenic or idiopathic [13]. These latter categories 
are here quickly reviewed. 

Clinically, pain can be labeled “nociceptive” if it directly correlates to the degree of 
tissue injury. More specifically, nociceptive pain is presumed to occur as a result of 
the normal activation of the nociceptive system by noxious stimuli. This type of pain 
response also occur in presence of inflammation (inflammatory pain). 

Neuropathic pain refers instead to dysfunction of the peripheral or central nervous 
system (CNS). These may be caused by injury to either neural or non-neural tissues. 
Although neuropathic pain is certainly influenced by ongoing tissue injury, it can be 
sustained by fundamental mechanisms which have become independent of the initial 
injury or damage. A continuing sensation of pain is in fact reported after the stimulation 
has ceased. Neuropathic pain presents diverse characteristics: On the one hand it mim- 
ics the qualitative features of somatic pain. On the other, it is also frequently described 
as a continuous burning, shock-like or paresthetic [14]. 

The patients psychological state contributes importantly to pain perception and as- 
sociated suffering symptoms [15]. The patient’ s self-report on pain should be evalu- 
ated while assessing other factors of paramount importance, as the presence of anxiety, 
depression, or other psychiatric disorders. In some cases, pain itself can be intimately 
linked to the psychological dimension, a phenomenon generically referred to as to “psy- 
chogenic” pain. However, when convincing inferences about the pathophysiological 
origin of the pain syndrome cannot be made, it is customary to diagnosis the so called 
“idiopathic” pain. 


12 Reducing pain via pharmacological therapy: Analgesics 


An analgesic, also known as a painkiller, is any member of the rather complex family 
of drugs, commonly employed to relieve pain or, equivalently, achieve analgesia. The 
word analgesic comes from Greek an- (“without”) and -algia (“pain”). 

Analgesic drugs act in various ways on the peripheral and central nervous systems. 
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The main classes of analgesics are the narcotics and the nonsteroidal anti-inflammatory 
drugs (NSAIDs).Th e narcotic analgesics, also termed opioids, are all derived from 
opium.Th e narcotic analgesics may vary in potency, but they are all extremely effec- 
tive when used in adequate doses. A number of chemical classes can be identified, but 
all share similar therapeutic effects and side effects. Notably, analgesics do provide a 
symptomatic relief, but, in general, they have no direct effect on the generating cause. 

Although pain syndromes may be different in many respect, the sensory pathway is 
ultimately solicited by the affected organ/tissue to the brain. Analgesics act at the level 
of the nerves, either by blocking the signal from the peripheral nervous system, or by 
distorting the subsequent processing by the central nervous system. 

To clarify the key elements that are to be encapsulated in a sensible mathematical 
model, as developed in the following, we shall here concentrate on reviewing the prin- 
ciple of opioids’ action. 


1.2.1 Opioid analgesics 


Beneficial and/or adverse effects of opioid analgesics can be traced to their interac- 
tion with the endogenous opioid systems [16]. Opioid agents and their receptors pop- 
ulate the central and peripheral nervous systems and other tissues. Opioid systems are 
indeed crucial in a vast gallery of homeostatic functions and movement control, as well 
as involved in the processing of noxious sensory input. The antinociceptive system, 
from which pain modulation stems, is itself extraordinarily complex. Information on 
this system constitute a useful background for understanding the effects associated to 
opioid analgesics. 

Pain transmission in the spinal cord is regulated by a balance of facilitatory and in- 
hibitory influences. These latter operate on the neural circuits of the somatosensory sys- 
tem. Noxious stimuli activate high-threshold primary sensory neurons in the periph- 
ery. This seeding activity is then conducted to their central terminals, which synapse on 
nociceptive (secondary) neurons in the spinal cord. Although opioid compounds are 
found in the periphery as well, they produce the condition for analgesia by primarily 
inhibiting the nociceptive transmission in the central nervous system (CNS). 

Opioid receptors positioned presynaptically and postsynaptically at the first cen- 
tral synapse in the spinal cord have been extensively analyzed. Those located on the 
presynaptic nerve terminal decrease the release of excitatory neurotransmitters from 
nociceptive neurons and respond to a large spectrum of noxious stimuli. Such a presy- 
naptic inhibition reflects by the opioid receptor activation on ion channels. More into 
details, opioid activation yields to hyperpolarization of the terminal through the open- 
ing of potassium channels or closing of calcium channels, the hyperpolarized neurons 
having a smaller probability to give rise to spontaneous discharge or evoked responses. 
Opioid receptors located postsynaptically have similar effects on the secondary neuron. 
Hyperpolarization caused by changes in ion fluxes leads to a reduced response of this 
neuron as an excitatory input is sent by primary nociceptive neurons. 

Opioids exert their analgesic effects by binding (and then activating) receptors that 
are part of the endogenous opioid system. This latter usually operates so to modulate the 
sensory input as caused by noxious stimuli, its response being activated by endogenous 
peptide neurotransmitters. Opioids can in turn mimic and, importantly, amplify the 
actions of those neurotransmitters. 

The endogenous opioid system includes a large plethora of opioid peptides which 
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are ligands for numerous types of opioid receptors. Some of these naturally produced 
peptides, induce a morphinelike effects and can be displaced from their apposite binding 
sites by opioid antagonists. Three distinct families of endogenous opioid peptides are 
often invoked: the endorphins, the enkephalins, and the dynorphins, which derive from 
the precursor polypeptides pro-opiomelanocortin, proenkephalin, and prodynorphin, 
respectively. 

The endogenous opioid peptides do bind to opioid receptors. In the CNS, there are3 
primary types of opioid receptors that mediate analgesia and are respectively labeled u, 
k, and ô. Enkephalins interact preferentially with the ô receptor, dynorphins with the 
k receptor, while endorphins bind to both 4 and ô receptors with a similar degree of 
chemical affinity. As previously remarked, those peptides display rather diverse physio- 
logic functions, one of which is associated to the antinociception effect here discussed. 
In different systems and settings, they can function as neurotransmitters, neuromodu- 
lators or, even, neurohormones. 

Drugs targeted to opioid receptors are broadly divided in four classes: agonists, par- 
tial agonists, mixed agonist-antagonists, and antagonists. Receptor activation by an 
agonist launches the pharmacologic actions, whereas an antagonist sit on the receptor 
without inducing any sensible effects. The intrinsic activity of a drug is then regarded 
as a quantitative indicator to distinguish between the aforementioned classes, includ- 
ing the intermediate categories. Partial agonists may also have antagonistic properties, 
because they compete with pure agonists for occupancy of opioid receptor sites. The 
degree of competition is determined by their affinity score to the receptor. For exam- 
ple, Buprenorphine hydrochloride, an analgesic used for addiction therapy, is a par- 
tial agonist being characterized by a very high affinity for the ju receptor; it can hence 
chase for the receptor and so have antagonist properties. The opioid analgesics most 
commonly used in clinical practice bind selectively to the u receptor and are called u- 
agonists. Morphine is the prototypical example of u-agonist. Despite the similarities 
between morphine and other u-agonists agents, different drugs can result in a large 
variety of effects in the individual patient. Consider for instance a patient who is chron- 
ically exposed to a -agonist and assume that it is suddenly switched to another: Pain 
can often be controlled by administering doses of the second drug that are by far lower 
than estimated according to their relative potencies and both the pattern and severity 
of non-analgesic effects can be distinct. This observation, widely known as incomplete 
cross-tolerance, suggests that these u-agonists are not acting through identical recep- 
tors. 

As it should be clear from the above, a consistent molecular model for drug absorb- 
tion requires accommodating for several different effects, which simultaneously cooper- 
ate with positive or negative interference. In the forthcoming chapters we shall focus on 
a limited set of key mechanisms which, we believe, do play a role of paramount impor- 
tance. In particular, anticipating our developments, we shall describe the dynamics of 
the active molecular agents by including the effect of the stochastic search for the target 
receptor and the competition with other chemical species populating the bloodstream. 


13 Importance of personalized medicine 


In medicine, individual response is crucial. Two patients, exposed to the same treat- 
ment can, for instance, experience completely different effect. Severe, even life-threate- 
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ning, side effects can be reported in one case, while in the other only minor conse- 
quences are registered. Alternatively, the drug may shrink an illness in one person but 
not in another. 

The body is a complex machine, and so many factors participate in making the re- 
sponse to an external input hard to predict. One major source of distinction between 
individuals stems from the inherited variations in the individual genes. Even slight vari- 
ations can alter the reaction of the body. Pharmacogenomics is the science that studies, 
among other things, how individuals react to an administered medications. Pharma- 
cogenomics is sometimes described as “personalized” or “individualized” medicine be- 
cause it holds the promise to eventually devise specific drug treatments based on the 
individual genetic background. 

Pharmacogenomics is indeed a promisingfi eld. A handful of tests are nowadays 
made available that can detect some of the genetic variations and so help predicting 
how a patient is likely to respond to an imposed medical treatment. Just to clarify, it 
could be revealed via the genetic analysis that a scrutinized patient is carrying a genetic 
modification which makes the drug stay in the body longer than normal, so eventually 
causing serious side effects. On the contrary, she/he may have a variation that reduces 
the effectiveness of the medication: This is hence less potent than reported on average. 
Once a variation is identified, scientists might be able to match it up with a response to 
a particular medication and so develop a personalized approach to medicine. 

Often the sensitivity versus a drug is due to a very small genetic modification that 
yields to a decreased activity of a particular enzyme responsible for the biotrasforma- 
tion of that drug. More precisely, once introduced in the body, drugs undergo a pro- 
cess of metabolization which converts them to metabolites. These latter are more wa- 
ter soluble with respect to the original compound and thus can be more easily ex- 
creted. Metabolism can also convert prodrugs (pharmacological substance adminis- 
tered in an inactive form) into therapeutically active compounds, and it may even re- 
sult in the formation of toxic metabolites [17]. Pharmacologists classify pathways of 
drug metabolism as either phase I reactions (i.e. oxidation, reduction and hydrolysis) or 
phase II conjugation reactions. Among enzymes that catalyze phase I drug metabolism, 
the most important are the cytochrome P450 enzymes, a superfamily of microsomal 
drug-metabolizing enzymes [17]. One member of this family, cytochrome P4502D6 
(CYP2D6) is by far the most intensively studied and best understood example of phar- 
macogenetic variations in drug metabolism. The CYP2D6 gene? is highly polymorphic 
with over 70 known alleles identified at the CYP2D locus on chromosome 22q13. At least 
15 of these alleles encode nonfunctional gene products as a result of single nucleotide 
polymorphisms (SNPs), gene deletion, aberrant splicing or premature translation ter- 
mination. Carriers of two nonfunctional alleles show a severely impaired metabolism 
of CYP2D6 substrates and are customarily referred to as to poor metabolizers (PMs). In 
contrast, individuals with at least one functional allele and thus normal CYP2D6 activity 
are termed extensive metabolizers (EMs). Among Caucasians, 5-10 % are PMs and fur- 
ther 10-15 % show impaired yet residual activity of CYP2D6, the so called intermediate 
metabolizers (IMs). The 1-5 % of the Caucasian population has a duplication or multi- 


> Genes encoding CYP enzymes, and the enzymes themselves, are designated with the abbreviation “CYP”, 
followed by an Arabic numeral indicating the gene family, a capital letter indicating the subfamily, and 
another numeral for the individual gene. The convention is to italicise the name when referring to the 
gene. 
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CYP2D6*1 Wild-type Normal 
CYP2D6*2 Several substitutions Normal 18 20 10 
CYP2D6*3 A2549 deletion Deficient 2 0 0 
CYP2D6*4 Gig4g A substitution Deficient 12-22 1-2 0-1 
CYP2D6*5 Gene deletion Deficient 2-7 4-6 6 
CYP2D6*9 G2613-A2615 deletion Decreased 2 2 3 
CYP2D6*10  CiooT substitution Decreased 12 4-6 51 
CYP2D6*7 C1023T, C2850 T substitutions Decreased 0 17-35 0 
CYP2D6*41  C_15g4G, G2988 A substitutions Decreased 8 10 3 
CYP2D6 x2 Gene duplication Increased 1-10 2-29 0-2 
CYP2C9*1 Wild-type Normal 

CYP2C9*2 C430T substitution Decreased 8-13 4 0 
CYP2C9*3 A1075C substitution Decreased 6-9 2 2-3 
CYP2C19*1 Wild-type Normal 

CYP2C19*2  GesgıA substitution Decreased 13 13-25 23-32 
CYP2C19*3  GesgA substitution Decreased 0 0-2 6-10 


Table 1.1: CYP allele subgroups, characteristic mutations, enzyme activity and frequency among Cau- 
casians, Africans and Orientals. Data are derived from references [18, 19, 20, 21, 22, 23] 


duplication of the CYP2D6 gene leading to the phenotype of ultra rapid metabolisers 
(UMs) [24]. 

Table 1.1 clearly indicates that the polymorphisms are not occasional mutations. On 
the contrary, they are widespread in the population. It is hence natural to assume that 
this pronounced degree of variability might interfere non trivially with the process of 
drug absorption. Before proceeding with a pharmacological treatment, one should have 
a clear picture on the genetic characteristic of the patient being treated. 

In the next section we will make reference to the case of a synthetic opioid, the tra- 
madol, to illustrate how the genetic polymorphism of cytochrome P450 deeply affects 
the mechanism of the drug’s action. 


1.4 The case of Tramadol 


In this section we review the main features of tramadol hydrochloride (tramadol), 
a centrally acting analgesic that is structurally related to codeine and morphine. The 
interested reader may refer to [25] for an exhaustive account on the subject. 

Tramadol was first synthesized in 1962 and has been available for pain treatment in 
Germany since 1977. It is administered as drops, capsules and sustained-release for- 
mulations for oral use, suppositories for rectal use and solution for intramuscular, in- 
travenous and subcutaneous injection. After oral administration, tramadol is rapidly 
and almost completely absorbed. It is mainly excreted via the kidneys. Plasma protein 
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Figure 1.1: Some known metabolic pathways of tramadol. 


binding is about 20%. 


The human in vivo metabolism of tramadol is complex with 23 metabolites be- 
ing identified: 11 phase I metabolites and 12 phase II conjugates. The major metabolic 
pathways (see Fig.1.1 ) are O-demethylation* to O-desmethyl-tramadol (M1) by the 
polymorphic isozyme cytochrome P450 2D6 (CYP2D6), and N-demethylation to N- 
desmethyl-tramadol (M2) by cytochrome P450 2B6 (CYP2B6) and cytochrome P450 
3A4 (CYP3A4). The primary metabolites of tramadol, may be further metabolized to 
three secondary metabolites, namely N,N-desmethyl-tramadol (M3), N,N,O-tridesme- 
thyl-tramadol (M4) and N,O-desmethyl-tramadol (M5). 


As already mentioned in the preceding section, the gene encoding for cytochrome 
CYP2D6 is known to show polymorphisms and the existence of different alleles results 
in functionally different enzymes. For this reason, the biotransformation of tramadol 
changes within the phenotypic population depending on the genotype of CYP2D6. More- 


4 Demethylation is the chemical process resulting in the removal a methyl group (CH3) from a molecule. 
In biochemical systems, this process is often catalyzed by an enzyme such as one of the cytochrome P450 
(CYP) family of liver enzymes. N-demethylation and O-demethylation are reactions which remove a 
group CH3 from NCH3 and OCH3. 
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Drug Opioid receptor affinity Uptake inhibition 

H ô K norepinephrine serotonin 
(+)-Tramadol 2.1 57.6 42.7 0.78 0.9 
(+)-Tramadol 13 62.4 54.0 2.51 0.53 
(-)-Tramadol 24.8 213 53.5 0.43 2.35 
(+)-M1 0.0034 
Morphine 0.00034 0.092 0.57 inactive inactive 
Imipramine 3.7 12.7 18 0.0066 0.021 


Table 1.2: Relative activity for inhibition of opioid receptor binding or monoamine uptake. Data (ex- 
pressed in pmol/L) from [26, 27]. 


over tramadol is administered as a racemic mixture’ of two enantiomers’, (+)-tramadol 
and (—)-tramadol, that are essentially metabolized by the liver producing (+)-metaboli- 
tes and (-)-metabolites, respectively. 

In vitro and in vivo studies have convincingly shown that the metabolism and dis- 
tribution of tramadol are stereoselective’. In other words, out of two (or more) possible 
reactions, one predominates. In vitro, O- and N-demethylation of tramadol were both 
shown to be stereoselective. The O-demethylation of tramadol, leading to M1, was deter- 
mined to be two fold greater for the (-)- enantiomer than for the (+)-enantiomer. On 
the other hand, N-demethylation, leading to M2, was considerably faster after incuba- 
tion of the (+)-enantiomer compared with the (-)-enantiomer. Since O-demethylation 
is the preferred pathway for biotransformation of tramadol, higher plasma concentra- 
tions of (+)-tramadol and (-)-M1 compared with (—)-tramadol and (+)-MI, respec- 
tively, can be expected in vivo. 

To study the stereoselectivity of renal clearance, isolated kidneys of rats were per- 
fused with tramadol and M1. The renal clearance of the enantiomers of both compounds 
was stereoselective, (-)-tramadol and (+)-MI being preferentially eliminated. In addi- 
tion, the O-demethylation of tramadol was stereoselective in the kidneys, (-)-tramadol 
being preferentially metabolized. 

Both tramadol and metabolites contribute to the analgesic activity via different mech- 
anisms. Tramadol displays only a modest affinity for u opioid receptors and no affinity 
for ô or K opioid receptors. The affinity of tramadol for u opioid receptors is approxi- 
mately 10-fold less than that of codeine and 6000-fold less than that of morphine. These 
scores for the affinity, taken as such, do not seem sufficient to explain the observed anal- 
gesic action of tramadol (see Fig. 1.2). 


In chemistry, a racemic mixture, or racemate, is one that has equal amounts of left- and right-handed 
enantiomers of a chiral molecule. 


Enantiomer is one of two stereoisomers (molecules that have the same molecular formula and sequence 
of bonded atoms, but which differ in the three dimensional orientations of their atoms in space) that are 
non-superposable complete mirror images of each other. 


Stereoselectivity is the property of a chemical reaction in which a single reactant forms an unequal mixture 
of stereoisomers during the creation of a new stereocenter. The selectivity arises from differences in steric 
effects and electronic effects in the reactions leading to the two products. 
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On the contrary, the metabolite M1 binds with about 300-fold higher affinity than 
the parent compound. Still, it has a much lower affinity when compared to morphine. 
The increase in subjective (and objective) pain thresholds as induced by tramadol is, 
at variance with what happens for other opioids, only partially blocked by the opioid 
antagonist naloxone. Therefore, the activation of u opioid receptors appears to be only 
one of the components of the mechanism of action of tramadol. More precisely, (+)- 
tramadol has a 2-fold higher affinity for the ~ opioid receptor than the racemate. Of 
the metabolites, (+)-MI has the highest affinity for the u opioid receptor, being about 
700-fold more potent than the (+)-tramadol. Another metabolite with a higher affinity 
than the (+)-tramadol for the u opioid receptor is (+)—M5. 

In addiction to its opioid action, tramadol inhibits the neuronal re-uptake of nore- 
pinephrine. In particular both the enantiomers of tramadol are involved in this pathway, 
the (+)-enantiomer being about 4-fold more potent than the (-)-enantiomer. More- 
over (+)-tramadol and its (+)-enantiomer, but not the (-)-enantiomer and MI, in- 
crease serotonin efflux. 

In conclusion, we can summarize saying that (+)-MI acts as a u opioid agonist, 
(+)-tramadol inhibits serotonin re-uptake and (—)-tramadol inhibits norepinephrine 
re-uptake. The activity of the other kind of metabolites has not yet been studied. 


1.5 Experiments 


In the preceding sections we presented a general, though synthetic, overview on of 
the main phenomena that are involved in pain expression and management. Those as- 
pects can be brought into evidence via dedicated experiments which enables researches 
to gain a comprehensive and quantitative insight into the problem at hand. 

When dealing with pain, it is crucial to quantify its associated intensity. Clinicians 
dispose of several, carefully validated, pain scales. Among them, the Visual Analog Scale 
(VAS) is by far the most adopted in medical practice. VAS is essentially a visual ques- 
tionnaire. The patient has to mark a sign on a graduated line (usually ten centimeters 
long) in the position that he feels would best correlate to the sensation of pain that he is 
experiencing. The reference positions are respectively the leftmost extreme on the line, 
labeled with “no pain’, and its end point which is assumed to correspond to the “worst 
possible pain” (see Fig. 1.2). 

Clearly, such a method cannot be regarded as an objective criterion for pain assess- 
ment, actual responses being strongly influenced by diverse psychological factors, as the 
emotional state or the anxiety. To overcome this obvious limitation, other methods have 
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been developed and extensively employed to quantify the effectiveness of a pharmaco- 
logical treatment. One of such objective criterion, often referred to as a physiological 
marker of pain, is the Evoked Response Potential (ERP). ERP is a time dependent elec- 
trical potential which is recorded from a human patient (or an animal), following an 
external stimulus. The time series of the brain electric activity is acquired via electrodes 
positioned on the head. The signal refers to the local activity of simultaneously firing 
neurons, and hence return a convolved information which necessitate further elabo- 
ration. Working with animals, more invasive tests have been conducted, allowing to 
unambiguously relate the registered signal to an externally imposed input. Aiming at 
clarifying the above technique, we shall here briefly recall the experiments described 
in [28]. In this paper, the authors investigate the effect of common anesthetic, on the 
rodent whisker sensory system. To this end an electrode array is implanted on the so- 
matosensory cortex of adult rats. The electrode is uniquely sensitive to the cluster of 
neurons which are connected to a specific whisker barrel. In this respect, it allows to 
register EEG and ECG signals (and their time evolution) after a whisker stimulation 
event. Figure 2 of [28] reports on the main conclusion of the aforementioned the study. 
Different panels represent the EEG and ECG amplitude respectively under anesthetized, 
wake, and sleep conditions. It can be immediately remarked that, while ERPs from con- 
trol (non treated) rats present rather stochastic patterns, regular cycles are instead found 
to emerge when the animal are kept under pharmacological treatment. This observa- 
tion constitutes indeed an indirect signature of medicaments’ action, which ultimately 
stems from their associated molecular peculiarity. As already noticed above, anesthetics 
act on large sets of biochemical reactions, being in principle directed towards indepen- 
dent neuronal families, and so determining profound differences in the evoked response 
components. 

Besides estimating the degree of experienced pain with a certain level of accuracy, it 
is also important to access a genetic screening of the patient being treated so to optimize 
the administered pharmacological protocol. Human DNA sequencing opened up novel 
scenarios and, among other things, translated into reliable methodologies for character- 
izing gene expression levels and detecting the presence of SNPs or deletion/duplication 
in the genome. 

These methods are essentially based on the well known microarrays technology. Mi- 
croarrays are solid surfaces on which a series of thousands of microscopic spots of DNA 
oligonucleotides are arrayed. Each spot contains a small portion of DNA sequence spe- 
cific for a gene. Th is can be a short section of a gene or other DNA element that are used 
as probes to hybridize a cDNA or cRNA sample previously labeled with a fluorophore 
dye.Th e idea is to determine relative abundance of nucleic acid sequences in the target, 
by monitoring the intensity of the emitted fluorescence. 

This is a rather powerful technology which enables for massive (and parallel) studies 
of genes activation within a given biological system. In the last few years, microarrays 
have been widely applied to thefi eld of genetic polymorphism. As an example, the Am- 
pliChip CYP450 Test manufactured by Affymetrix for Roche Diagnostics made it pos- 
sible to carry out a comprehensive analysis of two, particularly important, genes which 
encode enzymes crucial in mediating drug efficacy and adverse drug reaction. ‘The tests 
allow one to detect genetic variations in the Cytochrome P450 2D6 and 2C19 genes and 
provides the associated predictive phenotype profile (poor, intermediate, extensive, or 
ultra-rapid metabolizer). The AmpliChip CYP450 Test distinguishes among 29 known 
polymorphisms in the CYP2D6 gene, including gene duplication and gene deletion. 


12 Francesca Di Patti 


With reference to the CYP2C19 gene, the test can resolve the two major polymorphisms 
so far isolated. The actual implementation of the methodology is quite straightforward 
and non invasive. Only a sample blood from the patient is required. The blood is then 
manipulated in the laboratory so to perform a PCR amplification of selected segments 
of the DNA extracted from the patient sample.Th e amplified DNA segment is labeled 
with a fluorescent dye and then inserted into the AmpliChip CYP450, a plastic cartridge 
housing the microarray. The following step consists in the hybridization and scanning 
of the chip. The hybridization takes place into the Affymetrix hybridization chamber: 
The complementary base-pairs from the DNA fragments in the sample hybridize with 
those on the microarray which give a perfect match. Technical problems related on this 
specific operation will be reviewed in the forthcoming chapter. Then the Amplichip 
is moved to the Affymetrix scanner, where laser scanning of the hybridization pattern 
is performed. Images captured with the scanner are analyzed through a dedicated soft- 
ware which brings into evidence the genetic variations and quantifies the corresponding 
phenotype. 

The phenotypic profile of the patient constitutes a useful indicator which can help 
physicians to design the most appropriate drug and dose selection. Moreover, beyond 
the realm of purely genetic investigations, pharmacokinetics can also guide on the choice 
of the optimal therapy. The principles of pharmacokinetics, synergistically integrated 
with a pharmacodynamics viewpoint, are in fact useful ingredients to elucidate the com- 
plex dose-effect relationships. Figurel .3 presents a schematic interpretation of the in- 
terplaced complementarity referenced above. Pharmacodynamics correlates the con- 
centration of drug as measured in the blood to its induced effect. Pharmacokinetics 
aims at quantifying, under different conditions, the administered dose which is neces- 
sary to produce the sought concentration amount (as required by pharmacodynamics 
paradigms) [29]. This latter task necessitates accounting for various important ingredi- 
ents as absorption, distribution, metabolism and excretion (ADME). Absorption is the 
process through which the substance is propagates from the administration site to the 
systemic circulation. Distribution studies instead the dispersion of the drug molecules 
throughout the fluids and tissues of the body. Metabolism relates to the biotransforma- 
tion of the substance which yields to the metabolites. Finally, excretion clarifies the end 
process of the drug elimination from the body. 

Quantifying the four levels of the ADME process, implies estimating the associated 
pharmacokinetic reaction parameters. To this end a suitable mathematical model can 
be numerically fitted to the plasma drug concentration data as seen from blood samples. 
Textbooks on pharmacology [30, 31] reports indeed on very crude mathematical formu- 
lations which are by far reminiscent of the key physiological aspects involved. Indeed, 
rough compartmental ansatz are put forward which allow to define averaged mathemat- 
ical indicators. This latter are then extracted from real data via a simplistic inspection 
on the acquired experimental series. No serious attempt to perform direct, non lin- 
ear, fitting is registered that could provide a more sensible ensemble of the estimated 
parameters. 
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Figure 1.3: The relationship between dose and effect can be divided into pharmacokinetic (dose- 
concentration) and pharmacodynamic (concentration-effect) components. The concentration pro- 
vides a link between pharmacokinetic and pharmacodynamic domains: It is the target quantity for 
standard dose assignments. The three primary processes of pharmacokinetics are absorption, distri- 
bution, and elimination. 


Chapter 2 
Characterizing the individual response to medical treatment 


This chapter is devoted to presenting a selection of methods that we have developed 
to analyze different types of experimental data, as reported in the relevant medical lit- 
erature. We shall here discuss their potential interest when aiming at quantifying the 
individual response to a specific medical treatment. 

In the first section a deterministic model for the kinetics of tramadol is introduced. 
The model represents an example of a focused modelization grounded on the specific 
mechanisms which are essential elements of the dynamics of tramadol. As we shall see, 
by adjusting the theoretical curves to the experimentally measured concentrations, one 
can access an estimate of the chemical parameters which are then critically evaluated 
with reference to the metabolic profiles of the patients. The final aim of our work is 
to eventually obtain a quantitative indication on the required initial dose amount as a 
function of the phenotypic profile associated to the polymorphisms of CYP2D6. 

Our deterministic formulation moves from the observation that most of the known 
polymorphisms of cytochrome CYP2D6 are equipped by a phenotypic metabolic pro- 
file. Characterization of other isoforms of cytochrome CYP450, for example the isoform 
CYP3A4, which metabolizes tramadol, have not yet completely investigated, and only 
some hypothetical enzyme activity has been reported’. It should be emphasized that 
scientists have been interested in characterizing the cytochrome CYP2D6 because this 
latter is involved in the metabolization of nearly every drugs. At variance, cytochrome 
CYP3A4, appears to be active along a few metabolic pathways, and, for these reasons, 
results less attractive. 

Waiting for a detailed characterization of the cytochrome, one could, in principle, 
try to clusterize the genetic polymorphisms in order to hypothesize similarities in the 
corresponding enzyme activity. Standard clustering algorithm may result too rigid for 
being applied to this context, where, instead, one would need to dispose of efficient 
algorithms, suitable for the complex and variable concept of phenotypic distinction. 
Inspired to this rationale, we have developed an innovative procedure which is able to 
perform such tasks while resolving intermediate levels of the cluster structures. Section 
2.2 presents a detailed account on our implementation. Despite the method is at present 
not optimized for sorting real data on genetic polymorphisms of human cytochromes, 
we do believe that it could represent a rather promising starting ground for defining a 
future platform of data screening. 


1 The reader can compare the updated news relative to CYP2D6 and CYP3A4 at the Human Cytochrome 
P450 (CYP) Allele home page (http:www.cypalleles.ki.se). 
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The last section of this chapter copes with a common problem of the microarray 
experiments, namely the missing values problem. The new method that we propose is 
inspired to analogous developments in the broad field of opinion formation. Shortly, 
it consists in measuring the distance among records based on the correlations of data 
stored in the corresponding database. 


2.1 A comprehensive dynamical approach: Modeling the pharmacokinetics of 
Tramadol 


As introduced in the first chapter, one of the new challenges of the post-genomic 
era is represented by personalized medicine.Th is concept arises from the evidence that 
the same administered dose of drug can result in a wide variability of effects. These 
differences can be caused by acquired or inherited variability of absorption, distribution, 
metabolism and excretion of a drug. Pharmacogenetics and pharmacogenomics are 
emergingfi elds which aim at processing these information so to develop innovative 
medical protocols targeted to the individual patient. Pharmacogenetics is in particular 
concerned with detecting the genetical variability in the metabolism of drugs, which 
might yield to adverse drug reaction, toxicity or therapeutic failure of pharmacotherapy. 
Pharmacogenomics techniques are instead meant to isolate new pharmacological agents 
on the basis of a detailed knowledge of the human genome [32, 33]. 

As emphasized in section 1.3, any attempt to a personalized pharmacological treat- 
ment should account for the pharmacogenetics of the drug metabolism. Many drugs, in 
fact, undergo a metabolic transformation that produces active metabolites, occasionally 
more effective than the parent drug. Numerous genetic polymorphisms in drug metab- 
olizing hepatic enzymes have been reported and thoroughly characterized, in particular 
those belonging to the cytochrome P450 (CYP) superfamily, like cytochrome P450 2D6 
(CYP2D6). 

Moreover, we have already stressed that the conventional studies on the pharma- 
cokinetics of drugs are often restricted to analyzing the recorded data set on drug con- 
centration, bringing into evidence possible connection with pharmacodynamic aspects, 
such as the clinical response. These studies are in general not based on a rigorous math- 
ematical modelization, fully justified fromfi rst principles. On the contrary, they tend to 
involve rather approximate formulae which are then benchmarked to the macroscopic 
traits of the concentration changes over time (e.g. position of the peaks). As we shall 
clarify in the forthcoming discussion, we here take a different viewpoint by identify- 
ing a selection of relevant microscopic processes that constitute the backbone of our 
formulation. In doing so we aim at resolving, at least partially, the intricate dynamical 
interplay which drives the evolution of the interacting (chemical) species. This would 
in turn allow us to match the experimental data via a fitting procedure and so returning 
a quantitative estimate of the main parameters entering the model. These latter are then 
inspected as function of the, independently assessed, patients genetic variability. 

Before turning to introduce the model, it is useful to recall again the characteris- 
tic of tramadol. Tramadol (hereafter T) is a synthetic opioid, commonly used in the 
treatment of acute and chronic pain, which is known to be metabolized by CYP2D6. 
It acts as central analgesic and follows a complex pathway which is not yet fully eluci- 
dated [25]. However, it is generally believed that tramadol has a weak affinity with the 
4 opioid receptors, and it works through modulation of the noradrenergic and sero- 
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tonergic systems. Tramadol is administered as a racemic mixture of two enantiomers, 
(+)-T and (—)-T, which undergo hepatic metabolism via the cytochromes P450 form- 
ing five known metabolites [34]. Among the active ones, O-demethyltramadol (M1) is 
the most significant since it has 200 times the u-affinity of (+)-tramadol. The isoen- 
zyme responsible for O-demethylation is CYP2D6, that leads to the formation of (+)- 
O-demethyltramadol ((+)-M1) and (-)-O-demethyltramadol ((—)-M1). 


This cascade of reactions results in an underlying network of interactions which is 
still object of investigations. Once administered, tramadol diffuses in the blood cir- 
culation, where about 20% of the drug binds to the plasma proteins. Then tramadol 
molecules reach the liver where a fraction of it is metabolized by CYP2D6. The remain- 
ing amount abandons the liver and follows the circulation flux.Th e produced metabo- 
lites can be released in the circulation or, alternatively, may be further metabolized to 
secondary metabolites. Bloodstream drives the parent drug and metabolites to reach 
the target sites where they interact with receptors to produce analgesia. Elimination 
from the body is mainly due to the kidneys. As we shall see, the mathematical model 
here developed allows to track the time evolution of the main concentrations. 


It is worth mentioning that the applicability of the proposed model extends beyond 
the case under scrutiny. In fact we here focus on rather general physical (biological) 
mechanisms that are to be considered as key elements in any grounded pharmacokinet- 
ics theory. In this respect our approach is fully predictive and the observed evolution of 
the concentration amount are reproduced on the basis of selected interaction patterns. 
This is at variance with previously proposed scenarios [35, 36] where analytical laws for 
the concentration time curves are a priori guessed. 


We will proceed as follow. In section 2.1.1 we provide a short account on 1D spatial 
diffusion, a fundamental mechanism that certainly drives the process of homogeniza- 
tion and distribution of the drug in the blood. In section 2.1.2 we put forward our model 
for tramadol administered through an intravenous injection. The model results in a pair 
of coupled differential equations that govern the self-consistent evolution of tramadol 
and metabolites. Section 2.1.3 is focused on the validation of the model. First we present 
the experimental techniques and the measured data and then the theoretical model is 
fitted versus the measurements. Finally, we draw our conclusions and discuss a quanti- 
tative strategy for personalized drug treatment based on ourfi ndings. 


2.11 A short account on one-dimensional diffusion 


Once the drug is administered, for instance via intravenous infusion, it is trans- 
ported by the bloodstream and eventually reaches the target tissue. This is a typical 
passive transport process, the random motion of constituents resulting in a net change 
of the associated concentration, on a macroscopic level. Diffusion rules the dynamics 
and thus needs to be properly incorporated into the model. We will here review some 
fundamental aspects of diffusion theory that we shall further elaborate in the following 
section to derive a phenomenological entry to the problem of drug transport. 


Consider the case of a compound that diffuses along the x coordinate and suppose 
the distribution in the y — z plan to be uniform. We further assume a one-dimensional 
drifting velocity v in the x direction. The concentration at time t and position x, here- 
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Figure 2.1: The concentration p( L, t) is plotted versus time. Different curves refers to different param- 
eters (see legend). Here Sy, = 1. 


after termed p(x, t), obeys to the following transport equation: 
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where D stands for the diffusion coefficient.Th e previous equation can be analyti- 
cally solved, provided one specifies the initial and boundary conditions. For the sake 
of simplicity we shall here assume that initially (t = 0) the diffusing species is spa- 
tially located at x = 0, which, in our setting, ideally corresponds to the injection point. 
Mathematically, this amounts to require p(x,0) = d(x), where 6(-) stands for the so- 
called Dirac distribution. Notice that we have implicitly assumed the concentration to 
be normalized to one. As concerns the boundary conditions, it sounds reasonable to 
require that the concentration vanishes asymptotically, which formally translates into 
p(-oo,t) = p(co, t) = 0. Under these assumptions, the solution of equation (2.1) reads 


[37]: i 
7 1  (x-vt) 
Med = ER exp eno (2.2) 


where the cross-sectional area S,. of the system in the y — z plan has been introduced 
to convert the mathematical solution into real space (i.e. divide by the neglected di- 
mension). Focus now ona specific observation point, located at x = L, and monitor the 
time evolution of the concentration p(L,t). Equation (2.2) leads to: 


1 S (E) 
Sy2V4nDt i 4Dt | 


In Fig. 2.1 p(L, t) is represented versus time, for different choices of the parameters. 
As expected, the initial concentration is zero. A monotonic growth in subsequently 
displayed, and corresponds to the intuitive picture that an increasing amount of material 
finds progressively its way through the check-point at x = L. A clear peak is observed 


(2.1) 
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when the bulk material reaches the assigned spatial coordinates. The location of the 
maximum clearly relates to the diffusion properties of the medium (D,v) and to the 
choice of L. For larger times, the concentration diminishes and eventually fade out 
asymptotically, as prescribed above. 


Interestingly, the profiles reported in Fig. 2.1 bear an intriguing degree of similar- 
ity with the typical curves for the time evolution of drug concentration, as reported in 
a large number of experimental pharmacokinetics studies. This observation reinforces 
our belief that diffusion is indeed the key mechanism driving the process of drug ab- 
sorption. In a typical experiment, in fact, the substance is injected at a given location 
and then visits the intricate network of veins. Blood samples are then drawn from a 
distinct extraction point, and the operation is repeated at different times. By examin- 
ing the sample, one hence effectively monitors the time change of the concentration, 
as seen by a fixed reference control, a setting that closely resembles the scenario dis- 
cussed above’. In reality, blood diffusion occurs in a closed, topologically complex, 
loop. However, the tramadol data here analyzed do not present any clear evidence of 
the repeated (peaked) structures that one would expect in presence of a periodic flux 
[37], thus suggesting that the relevant dynamical process takes place in the time span 
that encompasses thefi rst passage.’ To simplify the discussion we shall therefore invoke 
a simple one-dimensional formulation, directly inspired to equation (2.1). Moreover, as 
we shall clarify in the following, the diffusion enters the proposed model as an external, 
time dependent, contribution. This is a minimalist approach that has the merit of en- 
abling us to formulate the dynamics in term of ordinary differential equations, without 
involving partial derivatives that need to be introduced in the framework of a rigorous 
description. 


2.1.2 The mathematical model: Diffusion and time delay 


As previously mentioned, the network of tramadol interaction is highly complex and 
a comprehensive picture that fully accounts for the underlying microscopic processes 
is still lacking. Here we isolate the main mechanisms in which T is involved, namely, 
the process of intravenous infusion, the biotransformation in M* and elimination, and 
build up a minimalist formulation which is shown to accurately interpolate the experi- 
mental data. A compartmental model is hypothesized and depicted in Fig. 2.2a.Th ese 
assumptions translate into the following differential equation for the concentration of 


Notice that also in the setting where injection and extraction points coincides, the drug has to ideally 
complete a full circulation tour, following the blood stream, before it can be detected in the sample drawn. 


Indeed, the paper by [37] cited above, solves the diffusion equation, in presence of drift, on a circular 
geometry. Within this setting, the time evolution of the resulting concentration, as measured at a given 
distance L, displays multiple peaks, which are originated from successive passages of the drifting pulse. 
For the case at hand, a direct inspection of the experimental data seems to exclude the existence of recur- 
rent patterns of the type mentioned above. In other words, it can be reasonably assumed that the drug gets 
dispersed, before the first complete loop is eventually closed up. Under these conditions, there is hence no 
practical difference in assuming the process to occur on straight line, as we do here, instead of including 
the details of the relevant periodic geometry. 


Here M incorporates all types of metabolites related to cytochrome CP2D6, including the M1 species for 
which data are available. 
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Figure 2.2: Schematic layout of tramadol (a) and metabolite M (b) dynamics.Th e local concentra- 
tion of tramadol in the blood, here labeled with T(t), changes as a result of diffusion through the 
circulation network.Th e time-dependent term F(t) is here introduced to account for such diffusive 
dynamics.Th e tramadol amount decreases due either to metabolization into M species (the rate of 
transformation being here km) or to generic losses (quantified through the reaction rate kz). The 
latter include both renal elimination and binding to receptors.Th e metabolites concentration M (t) 
increases, as result of the biotransformation of tramadol. Analogously to the above, M molecules are 
expelled by the body or bind to the receptors, thus resulting in an effective loss here controlled by the 
reaction constant h g. 


tramadol T(t), at the extraction point: 


Sr =-kmT(t)- keT(t)+F(t) (2.4) 
The term ky, T denotes the fraction of tramadol metabolized by CYP2D6 in M, while 
keT groups tramadol losses due, for instance, to the renal elimination, to the bind- 
ing of the molecules to receptors and to the formation of the other kind of metabolites. 
F(t) is an externally imposed, time-dependent, contribution which effectively mimics 
the role of diffusion, following the lines of the preceding discussion. The complex net- 
work that constitutes the pathway of blood circulation is here mapped into an idealized 
one-dimensional circuit. Practically, this assumption corresponds to neglect the occur- 
rence of periodic cycles of the diffusing substance, an event that can be excluded upon 
inspection of the recorded data. Mathematically, 


dp(L,t) 
dt 


where p( L, t) is given by equation (2.3). To labels the administered j1-mol of tramadol, 
thus providing the correct normalization. The drawing are assumed to be instantaneous 
and therefore we avoid to introduce further corrections to account for the time needed 
to complete the operation. 


F(t)=Ty 


Characterizing the individual response to medical treatment 21 


Metabolites M are produced from T, as schematically depicted in Fig. 2.2b. More 
specifically, at position L, we expect to measure an amount of metabolites which repre- 
sents the net balance between two competing processes. On the one hand, metabolites 
originates from a fraction of chemically active tramadol molecules, which are not barely 
transported by the above, purely diffusive, dynamics. On the other hand, metabolites 
may bind receptors and get also eliminated via the kidney, both mechanisms determin- 
ing a reduction of the detectable concentration M (t). One can hence hypothesize: 


where the delay 7 is here introduced to account for the finite time needed in the chemical 
conversion from T to M species”. 
Let us now focus again on the explicit expression for p( L, t). We shall introduce the 


variables a = L/V 4D and 8 = v/V 4D and recast p(L,t) in the form: 


a exp(- 2) 
Vro Tt t 


where Vrot = Syz L represents the total volume associated to the circulation network, 
Syz labeling here an averaged estimate of the cross-section surface. In the following we 
shall fit the above model to a set of experimental data, by properly adjusting the free 
parameters involved. By choosing the appropriate values that enable for an accurate 
matching between theory and experiments, we shall obtain a quantitative measure of 
the kinetic constants and eventually correlate our findings with the typology of patients 
under scrutiny. 

Finally, it should be emphasized that the simple model here discussed does not ac- 
commodate for a degree of negative interference among different species, (+)-T, (—)-T, 
(+)-M and (-)-M, a dynamical effect which arises when competing for available bind- 
ing sites (receptors). For this reason, we shall proceed with an independent analysis of 
data relative to the (+)-enantiomers and then to the (-)-enantiomers. The estimated 
parameters will be labeled with the corresponding symbols preceded by the prefix (+) 
or (-) (e.g. (+)-km stands for the kinetic rate of metabolites production relative to the 
(+)-enantiomers). 


p(t) = To (2.6) 


2.13 Validation of the model 


To validate the consistency of the model, we used data from [38]. Here we shall 
briefly discuss the experiments: Th e reader can refer to the original publication for fur- 
ther details on the procedures. In particular, we limit our discussion to the intravenous 
injection setting (termed phase C in [38]). Sixteen healthy volunteers, eight poor me- 
tabolizers and eight extensive metabolizers, took part in the clinical trial. The EMs hada 
median age of 26 years (range 25-28) and a median weight of 78 kg (range 71-85 kg); the 
PMs hada median age of 25 years (range 22-31) and a median weight of 76 kg (range 65- 
87 kg). The CYP2D6 phenotype was determined with sparteine and tramadol as probes, 


5 It is worth mentioning that the time delay 7 is an essential ingredient, which proves crucial to establish a 
quantitative agreement between theory and experiments. 
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Figure 2.3: Plots of experimental points (symbols) and numerical solution of the model (lines) for an 
extensive metabolizer. Plot (a) and (b) are relative to (+)-enantiomers, while (c) and (d) refer (—)- 
enantiomers. 


while Taqman technology was used to genotype four selected SNPs: the three known in- 
activation mutations CYP2D6*3 [39, 40], CYP2D6*4 [41, 42], CYP2D6*6 [43, 44] and 
the low activity allele CYP2D6*9 [45, 46].Th e volunteers received 100 mg tramadol 
hydrochloride in solution as an intravenous injection. Blood samples were drawn from 
an intravenous cannula in a forearm vein 0, 1/4, 1/2, 1, 1 1/2, 2, 3, 4, 6, 8, 10, 24, 34 and 
48 h after administration. The plasma samples were analyzed by the validated HPLC 
method [47]. 


Note that the (+)-tramadol time series from 3 subjects classified as PM were ex- 
cluded from the analysis. In those cases in fact the metabolites concentration was always 
found to be zero. 


The system of coupled differential equations (2.4) and (2.5) is integrated numerically 
using the matlab function dde23. The latter is designed to solve systems of differential 
equations with constant delays. In our scheme, the free parameters are tuned so to inter- 
polate the experimental time series for, respectively, the tramadol and metabolite con- 
centrations. A recursive algorithm is hence implemented to minimize an error function, 
defined as the sum of the differences between the theoretical and experimental curves, 
weighted according to their peak concentration to ensure a correct balance between the 
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Poor metabolizers Extensive metabolizers 


mean SD mean SD 
min] 2.5.20" 152107 > 19:10" 45-10” 
VJ/min] 1.7-10*° 1.3-107 1:9+10*° 2.8.107 
1/vimin] 8.1-1074 5.4-1074 7.9-1074 9.0-1074 
1/Vmin] 1.8510 2.6210" 1.1-108 1,8+107 
+)-vro [m8] 287108 «10108 2.19 1.4-10* 


-kp [l/min] 9.9-104 2.1-107 1.9-10° 1.7-10°3 
-)-kg [l/min] 2.0-10°? 1.5.1075 2.4.10 2.8.10 
+)-ky [l/min] 3.7.1075 2.0-10°° 5.2.1074 2.8-107 


m [l/min] 3.5-10* 1.5.1074 87-10% 4.0-10~ 


g [l/min] 1.9.10 16-10% 1.4-10° 1.4-107 
g [l/min] 1.1-10% 9.0-10* 2.0-10°% 1.9-1073 


) 
) 
) 
-)-vror [m8]  2.2-10*} 4.5-10*°  2.4-104! 1.4-10* 
) 
) 
) 
) 
) 
) 


Table 2.1: Average and associated standard deviation of the main parameters as obtained from the nu- 
merical fitting procedure. Th e (+) and (-) labels that precede the parameter symbol recall that values 
are calculated either from the (+)-enantiomers or the (-)-enantiomers experimental data set. No- 
tice that the parameters a and ( bear no immediate interpretation, being merely rescaled quantities 
introduced to simplify the fitting scheme. 


two contributions®. The values of the parameters a, 3, Vro» km, he and kg are updated 
until convergence. The procedure is repeated for different values of the time delay’ 7 
and the final global error monitored: The value of 7 that results in the smallest deviation 
is selected, and the corresponding set of fitted parameters stored. For each patient, one 
can therefore access an estimate of the constants that control the process of kinetics of 
(+)-tramadol (resp. (-)-tramadol) and (+)-metabolites (resp. (-)-metabolites). This 
procedure results eventually in quantitative estimates for the fitted parameters. Average 
of the best fit values for the classes here considered are enclosed in Tab. 2.1 together with 
the associated error. 

In Fig. 2.3 we report a visual comparison between measured points (symbols) and 
theoretical profiles (solid lines). The curves refer to an extensive metabolizer and clearly 
confirm the adequacy of the proposed scenario. 

The rates of metabolization and elimination play an important role in the kinetic 
of tramadol. A statistical analysis over a large population of patients would certainly 
contribute to shed light onto the crucial differences among PMs and EMs categories, 
and possibly to reveal the relative importance of the microscopic processes involved. 
Despite the fact that the statistical ensemble here analyzed is limited, we are however in 
a position to draw some general conclusion. 

Figure 2.4 reports the km values as resulting from the above fitting procedure for 
respectively (a) the (+)-enantiomers experimental data, and (b) the (-)-enantiomers 


é The code is constructed in matlab and can be made available upon request. The matlab function fminsearch 
is used to perform an unconstrained nonlinear minimization (Nelder-Mead). 


7 The minimal change in 7 occurs in step 1074. The best fit value results 7 = 0.001 for all analyzed cases. 
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Figure 2.4: Values of km obtained from numerical fit. Panel (a) refers to the (+)-enantiomers exper- 
imental data ((+)-km), while panel (b) to the (-)-enantiomers ((-)-km). The wire-framed rectan- 
gular patches define the region spanned by the data and it is solely drawn as a guide for the eye. 


data. The parameters km are plotted vs an integer number that labels the individual. 


Figure 2.4a displays a clear clusterization tendency of km: The two groups of sub- 
jects segregate, the extensive metabolizers being associated to a sensibly larger value of 
km. More quantitatively, the rate of metabolization in EMs is about 10 times greater 
than the corresponding analogue in PMs. This result agrees with the intuitive picture 
that the main difference between PMs and EMs is the ability to metabolize the drug. 
As clearly testified in Fig. 2.4b, such difference is less pronounced for (—)-T: The data 
cover patches which are only approximately disjointed. This finding points to the fact 
that biotransformation of tramadol varies within the phenotypic population of exten- 
sive metabolizers depending on the genotype of CYP2D6[ 48]. In order to obtain the 
evidence of a more clustered distribution, one should probably enhance the statistics, 
ie. accessing a larger set of experimental data. 


Notice that the above scenario is compatible with the scores of the t-Test. The null 
hypothesis, i.e. the assumption that the means relative to the compared ensembles are 
identical, is in both cases rejected, but with different p-values: The value associated to 
the (+)-km data set (0.003), is lower than the one extracted from (—)-km (0.007). This 
in turn implies that in the former case poor and extensive metabolizers are more clearly 
trackable to independent portions of the parameter space. 


Importantly, having accessed a quantitative estimate of the metabolization rate k m 
might enable us to define a therapy protocol which is tuned on the genetic profile of 
the patient. According to [49] in fact (+)-M1 molecules display higher affinity for the ju 
opioid receptor, being about 700-fold more potent than the affinity measured for(+)- 
tramadol. In this respect, metabolites molecules, simply M in the present discussion, 
seem to result in more effective chasers of target receptors, when compared to their tra- 
madol analogue. It can be consequently argued that the metabolite concentration in 
the blood is to be maximized, so to enhance the chances of a persistent and sensible 
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screening of the available receptor sites. Following this line of reasoning, we set down 
to calculate the initial dose of tramadol that is required to produce the same amount 
of metabolites in PMs and EMs individuals. To this end, we now disregard the trans- 
port of the drug? and restrict ourselves to a linearized version of the proposed model. 
For relatively short times t (tens of minutes), the metabolite quota M scales approxi- 
mately as km Tot, a formal relation which holds for both EMs and PMs patients with, 
respectively, km = Kee and km = ce . Using the values reported in the Tab. 2.1’, one 
immediately deduces an estimate for the relative dose of tramadol administered to EM 
and PM subjects, which is required to produce an identical concentration of metabo- 
lites, at a generic time t. More specifically, labeling with TP (resp. TY) the initial 
dose prescribed to EMs (resp. PMs), yields: 


TY = 14TE™ (2.7) 


for (+)-enantiomers and 
DE (2.8) 


for (-)-enantiomers!°. This conclusion quantifies the importance of including infor- 
mation on the kinetics of the two enantiomers when planning for personalized drug 
therapy. It is however worth emphasizing that the above conclusions are reached on the 
basis of a rather limited set of case studies. Increasing the number of patients, so to en- 
hance current statistics, represents a crucial leap forward, when aiming at confirming 
the correctness of our findings !. 

Moreover, we can speculate on the difference between (+)-ky and (-)-ky. Act- 
ing as an index of metabolization, kwm can be associated to stereoselective property of 
O-demethylation in favour of (—)-tramadol. In in vitro experiments O-demethylation 
of tramadol was determined to be greater for the (-)-enantiomer than for the (+)- 
enantiomer [50]. Our analysis indicates that (—)-ky, is always greater than (+)-km 
in every subject. This result confirms therefore the above scenario, thus providing an 
indirect evidence for the stereoselective property of O-demethylation. 

The distribution of h p is rather sparse. The values of (+)—h p range from 3.1-1074 to 
3.2-1073, while (—)-h g ranges from 3.5-107* to 4.5-1073. An exception is detected for 
subject 12, whose values of (+)-h g and (-)-hg are close to zero. This may correspond to 
a slow process of elimination, or can be alternatively interpreted invoking a low affinity 
of the M molecules to the receptors. 

For what concerns kpg, we could not identify a significant difference between poor 
and extensive metabolizers: (+)-k g is found in the interval [7.3-1074, 3.6-1075], while 
(-)-kg scans [3.4-1075, 6.7-10-3 ]. For subjects 11, 13 and 16 we found instead (+)-ky ~ 


8 The latter is instead a crucial ingredient when aiming at extract an estimate of the kinetic constants in- 
volved via numerical fit of the experimental data. 


° As documented in Tab. 2.1, the (+)-enantiomers (resp. (—)-enantiomers) case corresponds to kee = 
(5.2 + 2.8) - 1074 (resp. kag = kM = (8.74 4) - 1074) and km = KI = (3.7 + 2) - 107° (resp. 
km = kẸM = (3.5 + 1.5) - 1074). 


10 In both case the error on the estimated ratio T? M /TEM is less then 5%. 


1 It should be again stressed that the above conclusion are reached by putting forward several simplify- 
ing assumptions (small time regime, disregard the role of diffusion) which could in principle heavily re- 
flect on the final outcome. The estimated dose falls however within a reasonable range as e.g reported in 
http://en.wikipedia.org/wiki/Tramadol. 
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0 and for subjects 3, 9, 11 and 16( —)-kp = 0. We here recall that the term kpg controls 
the fraction of tramadol that binds to receptors, the quota eventually eliminated and the 
amount which is metabolized into other types of metabolites. To gain a comprehensive 
understanding on the overall process, and interpret these findings, it would be crucial to 
access the nociception and look for positive correlation with the low affinity of tramadol 
with receptors. 

Furthermore, we notice from Tab. 2.1 that, as expected, the entries relative to a, 3 
and vr; are rather sparse, the latter being essentially related to diffusion and transport 
properties. 

As a possible extension of the current work, we plan to develop a more detailed 
scheme of reactions that includes the distinction between (+)-enantiomers and (-)- 
enantiomers. Moreover, we stress the importance of incorporating into the theory the 
competition of (+)-T, (—)-T, (+)-M, (—)-M for the receptors’ binding sites, so to repro- 
duce the difference in affinity that have been reported in experiments [49, 25]. Another 
interesting extension concerns constructing a realistic, compartment-based model:This 
latter could eventually allow for a proper inclusion of the relevant diffusive mechanism, 
here incorporated via an effective term. 

Finally, we emphasize that M1 is just one out 23 known metabolites [51]. It would 
be therefore important to enrich the analysis by disposing at least of data relative to 
the second active metabolite produced through an hepatic cytochrome, namely N- 
desmethyltramadol (M2). To this end, it would be particularly attractive to perform 
experiments aimed at registering the concentration of T, M1, M2, for patients whose 
genotype has been determined, and monitor, at the same time, the perception of pain. 

In conclusion, as outlined above, the phenotype plays an important role in determin- 
ing the optimal medical strategy. However a complete characterization of the enzymatic 
activity and its relation with the genetic polymorphisms are far to be fully clarified. More 
extensive investigations are needed to gain a comprehensive understanding of the whole 
process. In the long run one will eventually dispose of a large gallery of complemen- 
tary information. Processing this information clearly implies developing robust and 
reliable methods for their classification. Clusterization algorithms constitute a viable 
tool, which certainly deserve further analysis. In the following section we shall discuss 
a novel procedure that we have outlined to perform such a delicate operation. 


2.2 Unbiased tools for data processing: Emergence of homogeneous clusters 


Cluster analysis is used to classify a set of items into two or more mutually exclusive 
groups based on combinations of internal variables. The goal of cluster analysis is to 
organize items into groups in such a way that the degree of similarity is maximized for 
the items within a group and minimized between groups. 

The need of data clustering is common in many problem. In particular, we can con- 
sider an abstract system represented as an unweighted graph, composed of N nodes. 
The nodes are connected by a variable number of links. A set of nodes is considered a 
cluster if, in some sense, the nodes belonging to it are more connected among them rather 
then among the rest of the system. 

The previous statement is rather vague: it is in fact very rare that one can unam- 
biguously identify an isolated cluster. In general, one needs to impose a cut-off, and 
therefore the cluster structure has a certain degree of arbitrariness. If one knows with 
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a sufficient precision the model underlying data generation, and therefore the expected 
data distribution, then the problem of clustering reduces to the minimization of distance 
between the measured and expected distribution, using the cluster number, coordinates 
and size as tuning parameters. In general, however, this is not possible, and so generic 
methods can only give a range of possible clusters, according to the variation of a given 
parameter. When there are extended regions of the parameters for which the cluster 
structure is stable, one may say that, for a given “coarse graining’, there is clustering. 
Let us consider the case of an adjacency matrix, with blocks with a given probability 
of ones (links) along the diagonal, and a smaller probability outside. Let us call these 
probabilities pı and po, respectively. In the case pa = 0 and p; sufficiently high (above a 
percolation threshold), the only clusters are those linking the nodes in the block. How- 
ever, as soon as po > 0, there are links connecting the previous clusters. In the limit of 
infinite sizes of the blocks, one gets a clear distinction between blocks, but for finite sys- 
tems, and without the possibility of performing ensemble averaging, fluctuations may 
well cause a site to be more linked to a “wrong” cluster than to the “right” one. In any 
case, even in the most ideal scenario, there is a certain degree of arbitrariness in de- 
ciding that the blocks form clusters. The situation is obviously worse when there is a 
continuum variation of the density of links. 


Most of clustering methods rely on performing walks on the graph, marking visited 
nodes or links, and measuring graph-distances. This is a rather heavy task, especially for 
large graphs. More efficient algorithms are probably those not based on path counting. 
Among these, one of the most promising clustering method in bioinformatics is the 
Markov Cluster Algorithm (MCL). Although this method is rather effective and really 
simple to implement, it lacks a precise physical meaning. 


Starting from MCL, we developed a new method based on a biological background. 
More precisely, we are interested in how tofi nd community structures with patches of 
populations in competition. We choose to refer to this method as to Ecological Cluster- 
ing (EC). 

In the next section we briefly describe the Markov Cluster Algorithm, while in the 
last one we introduce our model and compare it with MCL. 


2.2.1 Markov cluster algorithm 


The Markov Cluster algorithm (MCL) is a method to find community structure de- 
veloped by Stijn van Dongen in his PhD thesis [52]. The algorithm essentially simulates 
a process of diffusion on a graph. It works alternating two operations called expan- 
sion and inflation, on a N x N transition Markov matrix T obtained by dividing each 
element C;,; of the adjacency matrix by the degree of node j. The element 7;; of the 
stochastic matrix represents the probability that a random walker goes from node i to 
node 7. In thefi rst step, the left stochastic matrix of the graph is squared, so that ev- 
ery element of the resulting matrix gives the probability that a random walker, starting 
from node i, reaches j in 2 steps. The second step consists in raising each single entry 
of the matrix 7 to some power a, where a is a real value. This operation, called infla- 
tion, increases the weights between pairs of vertices with large values of the diffusion 
flow, which therefore have a high probability to belongs to the same community. The 
final step is to divide each elements of every column by their sum, such that the sum of 
the elements of the column equals one and a new stochastic matrix is recovered. The 


28 Francesca Di Patti 


corresponding algorithm for MCL is the following: 


Ci; 
Ti;(0) = 5, Cu 
Tt+At) = TO? 
Ti;(t+At) = Mc) 
Da (Ths (t+ At)) 


After some iterations, the process converges to a stationary matrix whose elements are 
either zero or one. According with the value of a, the graph represented by this matrix 
may be disconnected, and in this case its connected components are the communities of 
the original graph. Practically, the clusters are the few lines with some nonzero elements. 
The final partition is clearly dependent of the parameter a used in the inflation step. 
Therefore the method may converge to several different partitions and there is no way 
to decide which are the most meaningful or representative. 


2.2.2 Ecological clustering 


Our idea is to develop an algorithm tofi nd community structure based on a biolog- 
ical approach. In particular we start from a very simple model for the selection process 
of two species labeled A and B. We hypothesize that the number of elements of each 
populations, denoted by x and y respectively, obey to 


= ax°- ox 
y = by°-oy (2.9) 


where a and b denote the fitness value of A and B, while @ is a density limitation that 
keeps the total population x + y constant. The exponent c is the parameter which quan- 
tifies the non-linearity of the growth. It is easy to show (see for example [53]) that for 
c < 1 there is a stable mixed equilibrium between A and B even if one population is 
characterized by a growth rate bigger than the other one. For c > 1 there is an unstable 
mixed equilibrium between A and B, and the case where only one population is present 
is stable.Th is means that if, for example, the space is populated only by species A, B 
cannot invade also if it has a higher growth parameter. 

We are interested in investigating what happens if we simulate a similar scenario on 
an undirected graph with N nodes, where every node has enough room for N different 
species, all with the same fitness. Denote by P;;(t) the probability that at time t in the 
i-th node there is the j-th specie. We initialize the graph populating each node i with 
only individuals of species i, namely P;;(0) =0 Vi + j and P;;(0) = 1. Without any 
interaction among individuals, we can hypothesize a simple rule: At each time step the 
algorithm is based on two processes, namely diffusion and selection. During the process 
of diffusion, in each node i a fraction p of individuals decides to emigrate, where p e R 
and 0 < p< 1. They choose a link l according to C;;/ XY, Ci, where C is the adjacency 
matrix, and move to that node. The other 1 — p fraction of individuals remains in the 
original node. Then a new selection phase takes place, since nodes that received more 
than NN individuals are pruned, and nodes with less than N individuals are repopulated 
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by reproduction. This mechanism translates into the following algorithm: 


P,;(¢+ At) DIGHE 


Therefore, in the EC method we have two parameters, exponent a and diffusion p. 
Time also could be considered a parameter, since the convergence to the asymptotic 
status is enough slow that several structures appear during relaxation. However, here 
we concentrate only in the asymptotic state. 

To test this model and to compare it with MCL, we followed the idea described in 
[54] and we applied the algorithms to a network with a known fixed community struc- 
ture. For this purpose, we generated a network with N = 128 nodes, split into four 
communities containing 32 nodes each. Pairs of nodes belonging to the same commu- 
nity are linked with probability Pin, while pairs belonging to different communities are 
joined with probability Pout. The value of Pout is taken so that the average number of 
links a node has to members of any other community, Zout» can be controlled. While 
Pout (and therefore Zout) is varied freely, the value of Pin is chosen to keep the total 
average node degree, k, constant, and set to 16. 

Instead of focussing only in correctly classifying the nodes, we preferred to put at- 
tention on the ability to detect different levels of clusterization. For this reason, we did 
not test the sensitivity of the two algorithms measuring the quantities suggested in [54], 
but we preferred to monitor the entropy. 

We defined the entropy S(t) as 


Ep zklnz 
so. Eee 


where z;,(t) is the global distribution of k-th specie, namely z(t) = X; Tx;/N for 
MCL and z;(t) = X; Pik/N for our model. 

The entropy S(t) of both methods converges to a stationary level that assumes val- 
ues between 0 and 1. The steady state of entropy is 0 when only one big community is 
detected, while it is 1 when every node is classified as a single cluster. 

Fig. 2.5 shows the stationary entropy as function of the exponent a, relative to the 
two model tested on networks with different Pout. As we can see from this Fig. the MCL 
algorithm seems to work in a standard way. We can observe that there are two critical 
values, @ and a» that identify the regions in which the method finds different number 
of communities. More precisely, MCL converges to a single cluster for 1 < a < ay, 
it founds four clusters for a between a; and ao, and it detects 128 communities for 
a > Qz. Even though a; and a depend on the characteristics of the network, usually 
a, = 1.5 and ag = 2. It is interesting to notice that in this case, the parameter a acts 
as the parameter c in model (2.9): The value @ = 1, like c = 1, is a critical value which 
separates the region which leads to a stable mixed equilibrium among the species (128 
clusters), to the one where other equilibrium configurations are allowed. 

Although some of the previous features may be found also in our model, the EC 
algorithm exhibits an important difference. In fact we can notice that the stationary 
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Figure 2.5: Steady state of entropy S as function of the parameter a. Cyan dots are points relative to 
our model, while the magenta lines refer to the van Dongen’s model. For the plot in panel (a) we used 
a network with pout = 0.03, in panel (b) Pout = 0.04, in panel (c) Pout = 0.05 and in the last panel 
Pout = 0.06. In all cases, we used p = 0.95. 


values of the entropy are not only three as for the MCL. Consider, for example, Fig. 
2.5b. The entropy relative to our model follows the behavior of MCL only for small 
values of a. As a increases the entropy does not reach the boundary value 1 but assumes 
mostly two values, namely 0.1 and 0.28. To visualize how entropy is related to the 
final structure of communities, we can have a look at Fig. 2.6. This figure reports the 
stationary state of the network and the communities found by the two algorithm, for the 
same parameters of Fig. 2.5b, but with different values of a. The blue dots are the final 
structure of the network, while the continuous vertical red segments and the continuous 
horizontal green segments identify the communities found by EC and van Dongen’s 
method respectively. We can see that while the MCL converges essentially to three kinds 
of structures (one big cluster, 128 single clusters, and four clusters), the EC algorithm is 
able to detect intermediate clusters. In conclusion we can say that the main feature of 
EC algorithm is the ability to recognize different levels of clusterization. These levels are 
found monitoring the behavior of the entropy. 


In many applications bioinformatics deals with another crucial problem, the so called 
missing data problem. ‘This is reviewed in the following where a possible approach to 
overcome its limitations are discussed. 
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Figure 2.6: Blue dots represent the stationary structure of the network, green segments are communi- 
ties detected by the van Dongen algorithm, while the red segments are the ones found with EC model. 
All plots refers to the same network generated with Pout = 0.04, while a varies. In panel (a) a = 1.9, 
in panel (b) a = 2.0, in panel (c) a = 2.1, in panel (d) a = 4.0, in panel (e) a = 7.4 and in the last 


panel a = 8.0. 
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2.3 The problem of missing data: Filling the gaps 


As introduced in section 1.5 a DNA microarray is a collection of microscopic DNA 
spots of probes, commonly complementary to some region of a gene, arrayed on a solid 
surface by covalent attachment to a chemical matrix. DNA arrays are commonly used 
for expression profiling, namely monitoring expression levels of thousands of genes si- 
multaneously, or for comparative genomic hybridization. Despite their wild use, gene 
expression microarray experiments, frequently generate data sets with multiple missing 
expression values. This may be due to a production flaw, a mistake during the sensible 
procedure of preparation, or the final image may be corrupted or may have an insufh- 
ciently resolution. Some researchers overcome the problem proceeding with the anal- 
ysis omitting missing values and data relative to suspected spots, like spots with dust 
particles, irregularities or other bad features. 

The presence of missing values implies loss of information, and the practice to ignore 
this would lead to an unbalanced experimental design. Moreover, many algorithms for 
gene expression analysis require a complete matrix of gene array values as input, and 
may lose effectiveness even with a few missing values. 

Methods for imputing missing data are needed, therefore, to minimize the effect 
of incomplete data sets on analyzes, and to increase the range of data sets to which 
these algorithms can be applied [55]. Moreover, comparison between a “forecasted” 
value based on correlations in the dataset, and the measured one, can be considered a 
consistency “check” of the dataset itself. 

Here we propose a new method for inferring missing data which takes the cue from 
the field of opinion formation. In opinion formation, one can assume that one’s opin- 
ion on a certain item is given by the characteristics of the item, weighted by individual 
“tastes” The tastes result from past experiences, but they do not change abruptly from 
time to time. In principle, tastes can be decomposed into independent “dimensions” It 
is rather difficult to identify such dimensions, as testified by the limited success of mar- 
ket campaigns. However, it can be shown [56] that exploiting the correlations among 
the expressed opinions, it is possible to deduce the distance between the tastes of two 
individuals. 

We start from the original method described in [56] to derive a new way for mea- 
suring the distance among records based on the correlations of data stored in the cor- 
responding database entries. In this work the opinions expressed over a set of topics 
originate a “knowledge network” among individuals, where two individuals are nearer 
the more similar their expressed opinions are. Assuming that individuals’ opinions are 
stored in a database, the authors show that it is possible to anticipate an opinion using the 
correlations in the database.Th is corresponds to approximating the overlap between the 
tastes of two individuals with the correlations of their expressed opinions. Here we ex- 
tend this model to nonlinear matching functions, inspired by biological problems such 
as microarray (probe-sample pairing). We investigate numerically the error between 
the correlation and the overlap matrix for eight sequences of reference with random 
probes. Results show that this method is particularly robust for detecting similarities in 
the presence of traslocations. 

Let us first illustrate the problem summarizing the main results reported in [56]. 
Consider a population of M individuals experiencing a set of N products. Assume that 
each product is characterized by an L-dimensional array a = (al), a®,..., al) of 
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features, while each individual has the corresponding list of L personal tastes on the 
same features b = (6°) b)... bŒ). The opinion of individual m on product n, 
denoted by 5, n, is defined proportional to the scalar product between bm and an: 
Sm,n = A(L) bm: an, where A(L) is a suitably chosen normalization factor. In general, 
A(L) should scale as L~* and depend on the ranges of a and b. 

In order to predict whether the person 7 will like or dislike a certain product a,,, 
assuming to know an, it is sufficient to obtain the individual tastes of that individual, i.e. 
the vector b;. The similarity between tastes of two individuals è and j is defined by the 
overlap ();; = b; - bj between the preferences b; and b,. 

One can build a knowledge network among people, using the vectors bm as nodes 
and the overlaps (2;; as edges. Maslov and Zhang [57] (MZ) assume that a fraction p of 
these overlaps are known. They show that there are two important thresholds for p in 
order to be able to reconstruct the missing information. The first one is a percolation 
threshold, reached when the fraction of edges p is greater than pı = 1/M — 1 where M 
is the number of people. This means that there must be at least one path between two 
randomly chosen nodes, in order to be able to predict the second node starting from 
the first one. 

Since vectors b,, lie in an L dimensional space, and a single link “kills” only one 
degree of freedom, a reliable prediction needs more than one path connecting two in- 
dividuals. Maslov and Zhang show that there is a “rigidity” threshold po, of the order 
of 2L/M, such that for p > ps the mutual orientation of vectors in the network is fixed, 
and the knowledge of the preferences of just one person is sufficient to reconstruct those 
of all the other individuals. 

In general one does not have access to individuals’ preferences, nor one knows the 
dimensionality L of this space. In order to address this problem, the authors define the 
correlation C;; between the opinions of agents 7 and j by 


N è pix gT. È a arate 
Ce = En-1(Sin Si) (Sjn 5;) l (210) 
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where 5; is the average of the opinion matrix S over column i. The elements C’,; can be 
conveniently stored in a M x M opinion correlation matrix C. 
One can compute an accurate opinion anticipation S,,,, of a true value Smn using 


this formula: 
k M 


Smn av C misin 2.11 
Mi (2.11) 


where k is a factor that in general depends on L and on the statistical properties of 
the hidden components. However, if the components of an and bm are independent 
random variables, k is independent of n and m, so it can be simply chosen in order to 
have Šmn defined over the same interval as Smn. 

For large values of N and M, the factor k can be identified with the number of 
components L, and obtain an estimate for the average prediction error 


E= ui X (Smn = A x yr YM+YN (2.12) 
MN a VMN 
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where 
y= MLN}. (2.13) 


Formula (2.12) implies that the predictive power of equation (2.11) grows with M N and 
diminishes with L. This fact is a consequence of the decay of the correlations among 
opinions with L, so that more amount of information is needed in order to perform a 
prediction as L grows.Th is condition can be compared with the “rigidity” threshold pə 
in the MZ analysis. 


2.31 Test case microarray inspired 


In order to investigate the introduction of nonlinearities in the function used to 
model the process of opinion formation, we considered the case of a microarray. 

As mentioned before, microarray experiments can suffer from the missing values, 
and this fact represents a problem for many data analysis methods, which require a com- 
plete data matrix. Although existing missing value imputation algorithms have shown 
good performance to deal with missing values, they also have their limitations. For 
example, some algorithms have good performance only when strong local correlation 
exists in data, while some provide the best estimate when data is dominated by global 
structure [58]. 

Here we modified the model described in the previous section to investigate the 
relationship between the correlation and the overlap between sequences. 

To do this we considered an alphabet of four symbols, namely A, T, G, C, correspond- 
ing to the four nucleotides that constitute the DNA. We used this alphabet to generate 
randomly M sequences of length L representing the probes of the microarray. Then 
we generated N samples of length W representing the sequences to be hybridized on 
the microarray. 

The correlation C;,; between sample è and sample j is defined by 


Epi (Mik — Mi) (Mix — my) 


Cy = i,j=1,...,N, (2.14) 


Vili (Mik — Mi)? Dy (Myx - M)? 


where mx is the maximum complementary match between sample è and probe k with- 
out gaps. 

The aim is to test the relationship between the correlation matrix C and the overlap 
matrix Q constructed using the following idea of similarity. We hypothesized to infer 
the similarity between sequences based on the number of subsequences of length L in 
common. For this reason we defined the overlap (2;; between sequence 7 and sequence 
j as the number of subsequences of length L that appear in the both sequences, divided 
by W - L + 1 for normalization. This matching function is nonlinear since the effect of 
a mismatch depends on its position in the subsequence. 

To test our hypothesis, we considered eight referential sequences: 


Seq. 0: This is the first reference sequence, completely random of length W. 


Seq. 1: Equal to sequence0 , except for a mutation in the middle (this mimics the Af- 
fimetrix central mismatch mechanism for measuring the level of random pairing). 


12 The probes in real microarray are discriminated generally carefully chosen in order to genes of interest. 
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Figure 2.7:Th e error as a function of the length W of the samples, averaged over 40 realizations, 
N = 10, L = 30, M = 500.Th e plots of sequences 0-2 and 0-3 refer to the right y-axis. One can 
observe that all errors diminish with W. 


Seq.2: Equal to sequence0 , but shifted of one basis. 
Seq. 3: Equal to sequence0 , with shift and central mutation. 
Seq. 4: First half of sequence 4 is equal to the second half of sequence0 , and vice versa. 


Seq. 5: First half of sequence 0 is equal to the second half of sequence0 , the rest is 
random. 


Seq. 6: Another reference sequence. 


Seq. 7: Sequences 6 and 7 contains the same “gene”, of length W/3, in different posi- 
tions. 


To check the validity of the model described in the previous section, we measured 
the error ¢;; for the pair of sequences è and j defined as the absolute value of the differ- 
ence between the correlation and the overlap, namely £;; = |C;; — (;;|. We performed 
various simulations and than we calculated the average of the error denoted by £. 

In Fig. 2.7 we plotted the error € vs W. One can see that for all the analyzed cases 
the error decreases, and this result agrees with those reported in [56] (the parameter W 
here corresponds to M in the opinion formation model). 

Figure 2.8 shows the behavior of £ with respect to M. The curves are approximately 
constant, showing that the error is independent of M. 

As one can see from Fig. 2.9, where we plotted the error vs L, € does not follow 
a monotonous trend, except for the pair of sequences0 -4 for which the value of £ is 
almost constant and next to zero, and for €9; which increases. For what concerns the 
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Figure 2.8:Th e error £ as a function of the number of probes M, averaged over 40 realizations, N = 10, 
L = 20, W = 150. One can observe that errors do not vary with M. 
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Figure 2.9:Th e error € as a function of the length L of the probes, averaged over 100 realizations, 
N = 10, W = 200, M = 1000.Th e plots of £02, €03 and €04 refer to the right y-axis. 
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values of €02, €03 and Eos, one can detect that the errors decrease until L ~ 10, because 
probes too short can hybridize in many positions without a high specificity. Then they 
oscillate until L = 35, and for larger L the errors increase. This last increase is due to 
the small coverage of the probes in the sequence space, since we kept the number of 
sequences M fixed while the sequence space grows as 4”. 

To sum up, the monitoring of the error for the eight sequences of reference, with 
respect to M, W, and L, allows us to assert that the error is low in all cases, decreasing 
when W increase, and independent of M. With respect to L we find that the model is 
more robust for traslocation. 

In conclusion we can say that the correlation matrix of our model can be used to 
estimate the distance between sequences. Moreover we point out that the same result 
can be found following the idea of “negative database”, namely using the subsequences 
of length L not in common between two sequences.Th is idea arises from the antibody- 
antigens domain which is commonly affected by the missing values problem. 

Antibodies are proteins that are used by the immune system to identify and neutral- 
ize foreign objects, such as bacteria and viruses. Classifying antibodies, based on the 
similarity of their binding to the antigens, is essential for progress in immunology and 
clinical medicine. A striking feature of the natural immune system is its use of negative 
detection in which “self” is represented (approximately) by the set of circulating lym- 
phocytes that fail to match self.Th is suggests the idea of a negative representation, in 
which a set of data elements is represented by its complement set. That is, all the ele- 
ments not in the original set are represented (a potentially huge number), and the data 
itself are not explicitly stored.Th is representation has interesting information-hiding 
properties when privacy is a concern and it has implications for intrusion detection. 
One of the example where this idea has been concretized is exactly the case of a nega- 
tive database [59]. 

In a negative database, the negative image of a set of data records is represented 
rather than the records themselves. Negative databases have the potential to help pre- 
vent inappropriate queries and inferences. Under this scenario, it is desirable that the 
database supports only the allowable queries while protecting the privacy of individ- 
ual records, say from inspection by an insider. A second goal involves distributed data, 
where one would like to determine privately the intersection of sets owned by different 
parties. For example, two or more entities might wish to determine which of a set of 
possible “items” (transactions) they have in common without reveling the totality of the 
contents of their database or its cardinality. 

Since in our model a datum is essentially stored as the set of matching items plus 
the set of non-matching ones, our results can be applied both to positive and negative 
representation of data. 


Chapter 3 
Role of fluctuations in the experienced pain perception 


In the first chapter we have reviewed the main features of pain, focussing on the 
cascade of successive reactions which are activated by the so called noxious stimuli: 
The peripheral terminals of primary sensory neurons launch the signal, which is then 
transmitted to the spinal and supraspinal nuclei and eventually yields to the activation 
of a matrix of cortical areas that are deputed to the conscious experience of pain. 

More specifically, the stimulus originating from a bodily harming menace can be di- 
rectly processed through transduction by the receptors located on the nerve terminals. 
Alternatively, an indirect pathway can take over through the activation of transient re- 
ceptor potentials on keratinocytes or the release of intermediate molecules which, in 
turn, act on sensory neurons receptors. In the following we shall assume the first sce- 
nario to hold, and, though certainly important, disregard other mechanisms that might 
be simultaneously in play. In other words, we simplistically imagine that pain receptors 
act as effective gates, channeling the route to the involved cortical circuits. 

Analgesic drugs relieve the pain by interfering with the peripheral and central ner- 
vous system. Drug molecules bind in fact their target receptors, and consequently in- 
hibit the pain perception. To grasp and visualize the essence of the process, one can 
hypothesize that the bound chemical element occludes the path, by impeding the signal 
transduction through the channel envisioned above. 

Analgesic are commonly used in basic research and clinical practice, but their inter- 
action with nociceptory and normal sensory processing remains to be fully unraveled. 
In section 1.5 of the first chapter we mentioned that analgesics are for instance known 
to modify the electrical recordings measured via evoked potentials (ERPs) responses 
[60], a powerful diagnostic tools employed to monitor and characterize a large variety 
of central nervous system disorders. ERPs are elicited by a specific stimulus applied to 
the e.g. pain receptors and consist in recording the induced electrical brain activity, as 
detected by localized electrodes placed on the surface of the head. Furthermore, ERPs 
are also useful in documenting objective response to pain [61, 62] and can thus prove 
fundamental to elucidate the molecular processes that control analgesic absorption and 
metabolization. 

As shown by Fig. 2 of [28] discussed in section 1.5, different analgesic agents have 
been shown to produce intriguingly distinct effects at the level of the ERPs. Recorded 
time series of the solicited electric activity display in fact remarkably different patterns, 
which are generically attributed to the chemical specificity of the analgesic compound. 
Qualitatively, large, regular, oscillations of the electric response manifest, latency and 
amplitude being peculiar traits, supposedly related to the molecular characteristic of 
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the administered drug. 

Furthermore, cycles in the perception of pain have been also reported which might 
be hypothetically driven by similar microscopic processes, the interaction between the 
anaesthetic molecules and their targets playing certainly a role of paramount impor- 
tance. Clearly, the individual experience of pain is also influenced by psychological and 
cultural factors, unfortunately difficult to deconvolve when aiming at resolving the ob- 
jective picture. 

The issue of developing a unique interpretative framework to account for the pres- 
ence of such oscillatory regimes has catalyzed vigorous discussions. The puzzle of their 
existence remains however to be fully understood. 

Current mathematical models approach the problem via deterministic paradigms, 
thus neglecting the crucial role which is certainly played by the noise, intrinsic to the 
phenomenon under scrutiny [63].Th ese aspects become particularly important when 
accounting for the presence of diverse chemical species, which populate the stream flow 
in a spatially diffusive environment. Different chemical entities may compete with the 
drug molecules and occupy the sites located in close vicinity of the receptors, thus ef- 
fectively hindering the binding event. Under specific conditions, such competition sus- 
tained by the stochastic component of the dynamics might result in large temporal os- 
cillations for the amount of bound receptors, a mechanism which could explain the 
emergence of macroscopic cycles for the sensation of pain in response to medicaments. 

In this chapter, we shall speculate on the above scenario by putting forward a net- 
work of chemical reactions and performing a system-size expansion through the cele- 
brated van Kampen theory [64]. This enables us to derive a set of linear equations for the 
fluctuations, with coefficients related to the steady-state concentrations predicted from 
the first-order theory (i.e. the deterministic rate equations). Solutions are identified 
for which the deterministic steady-state occurs via damped oscillations: the inclusion 
of second-order fluctuations leads then to the amplification of sustained oscillations. 
These conclusions are briefly discussed with reference to the existing medical literature. 


3.1 The chemical equations governing the microscopic process 


Within the simplified scenario depicted above, we shall model the chemical inter- 
action between a large, though finite, number of drug molecules (analgesics), hereby 
termed T, and free receptors Rp which represent their binding target. Following a suc- 
cessful binding event, a molecule of the species T disappears, leaving an empty case, 
hereafter labeled E. The population of bound receptor Rr is in turn increased by one 
unit. These assumptions formally translate into the compact chemical notation: 


Ret? = Rea kB (3.1) 


where a stands for the associated reaction rate. The inverse reaction corresponding 
to the spontaneous detachment of the bound component may occur! with a certain 


1 We here assume that the free T molecule is still chemically active and can thus potentially chase for un- 
screened targets. This working hypothesis can be relaxed leading to conclusions qualitatively similar to 
the ones highlighted below. 
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probability 6°, which motivates the introduction of the dual relation: 


Roepe Rp+T (3.2) 


The analgesic molecules T surf in a densely packed environment and certainly ex- 
perience collisions with several other microscopic chemical entities, H, which populate 
the streaming flow. Binary interactions between H and T elements, can occur in the 
close vicinity of the receptors (Rp) location, potentially disturbing and eventually inter- 
fering with the binding event. As a result of an hypothetical collision, the active species 
T can be ejected by the spatial layer immediately adjacent to the receptor, leaving be- 
hind an empty case E. Beyond this effect, which stems from purely steric interactions, 
one has to account for possible chemical transformations, which might occur when in- 
dividuals from the H and T species encounter: The active chaser T can lose its ability 
to bind the target? and it is thus mapped into an inactive H molecule. To incorporate 
these effects into the proposed description we postulate the following interaction rules, 
which are loosely inspired to the predator-prey competition mechanism: 


Fip > H+H (3.3) 


o 


H+T > H+E (3.4) 


The values of o and y characterize the effectiveness of the interaction, which is in turn 
sensitive to the choice of the compound T. The idealized cartoons of Fig. 3.1 are aimed 
at visualizing the above reaction schemes. 

To complete the model we introduce an effective migration, by requiring that the 
T and H molecules can enter (resp. leave) the region deputed to the interaction. The 
latter assumption yields to the following set of chemical relations: 


T & E (3.5) 
E & T (3.6) 
H È E (3.7) 
E <> H (3.8) 


The population, namely the ensemble of elements belonging to an homologous species 
X, will be labeled in the following with the symbol nx. Notice that the number of 
receptors N1 = nr, + nr, and the total amount of molecules (including the empties) 
Nə =nr+ny+Nz are conserved quantities. This observation enables us to reduce the 
complexity of the problem by setting: 


Nr, =Ni- nr, Np = No-np-Nny 


In principle it would be extremely useful to dispose of experimental estimates for the reaction rates, so to 
define a realistic range of variability for the free parameters in the model. The most reliable data concern 
the so-called (equilibrium) affinity constant for the case of e.g. the Tramadol. Depending on the target 
receptor (and on the specificity of the chaser’s molecule) the affinity constant is reported to vary of a large 
amount which scans two orders of magnitude (from fraction of unity to hundreds) [25, 49]. 


Note that this can happen both due to a mechanical stress or via chemical combination of the colliding 
species, see for instance [30] where the plasma protein binding is discussed. For a specific application 
relative to the case of tramadol refer to [65]. 
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Figure 3.1:Th e main reaction schemes are depicted.Th e squares stand for the inactive species, while 
the circles represent the drug molecules. The model is then complemented with a set of additional 
reactions (see equations (3.5)-(3.8)), which accounts for the possibility that T and H enter (resp. leave) 
the region deputed for the interaction. 


In the following we shall use the vectorial notation n = (nT, NR, NH) to help keeping 
the mathematical developments compact. 

If we denote by nx the number of elements belonging to an homologous species 
X, it is reasonable to assume that number of receptors Ni = nr, + MR, and the to- 
tal amount of molecules (including the empties) N2 = nr + ny + Nz are conserved 
quantities. The inherent complexity of the problem is hence reduced by setting: 


Nr, = Ni -Npp ng=Nə-nr-ny 


so that the system is completely specified by the vector n = (NT, NR; NH). 

We assume that the two populations are homogeneous, where all individuals have 
the same chance of interacting. In this case, the probability of an events happening 
only depends on the probabilities of choosing the various elements from the total cor- 
responding population and the rates at which they take place. So we can define tran- 
sition rates T(n'|n) from state n to a new state n’, where n’ is one of the admissi- 
ble state defined by the chemical equations (3.1)-(3.8). In our case, n’ can only be 
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(nrl, ne, Fl, ny), (nr-1, nrp ny+1), (nrl, NR, NH) Or (NT, NR, nyl), 
and so just seven transition probabilities differ from zero, namely: 


T -1 1 = na 
(nr » Rp + NH n) N Ni 
Na-nr-Nny NR 
T(np+1,ng,- ny) = 8 = 
( T R H ) Nə Ny 
NH NT 
T -1 1 = y=- 
(nr SMR NH + n) N>» Nə 
NH NT NT 
T -1, À = Sa ae Ol ae 
(nr nry,ny|n) N, No N; 
Nə-nr-n 
T(nr+1,nr,-,nun) = m 2 ~ a 
2 
NH 
T ’ ; -1 = da 
(nr Rr, NH n) 2 N, 
N, — = 
T(NT, NR NH + 1\n) = M - To e 
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With the above transition probabilities the master equation for the probability P(n, t) 
to observe the system in the state n at time t, may be written as 


L Pint) = T(alnr+1,nr,-1,n4)P(nr+1,nr, -1,ny,t) 


+T(n NT l,nry SR l,ng)P(nr = l,nry + l,ny,t) 
+T(n|nr+1,nr-;ng-1)P(npr+1,ng,,ng-1,t) 


+T(n|n7 +1,nR,,27)P(n7+1,NrR,, NH, t) 
+T(n|n7 -1,nR,,NH)P(n7-1,NR,, NA, t) 
+T(n\n7, NR, ny +1)P(n7,nR,, NH + 1,t) 


+T(n\n7,nR,, NH -1)P(nr,nR,, NH - 1,t) 
-[T(nr - 1, nr, +1,ng|n)+T(np+1,nr,-1,ny|n) 


+T (nr 53 1, NRr; H + 1|n) + T(nr == 1, nr,, nun) 
+T(nr+1, nry, nun) +T(n7, TR) NH- 1|n) 


ED tie NH + Un) |P(n, t) (3.9) 


Multiplying both sides of equation (3.9) by nx, summing over all states, and remem- 
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bering that by definition (n) = Xn nx P(n, t), we obtain 


oie ) = em vo) (TATE) _ 5 (“=| 
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In the limit Ni, N2 — oo, we can replace (nyny) with (nx)(ny), and (n7)/Na, 
(ner)/N,, (nz)/ N; become the deterministic variables dr, der and dy of the mean- 
field system 


Lor = -G@r(1 Prr) + BOR, (1- br - bn) 


-(1+6)dndr è. dr + (1 or — bn) (3.10) 
Lón = e[@gr(1-6r2)-Ada:(1-9r-$n)] G0 
= 
Lon = uer +Î(l- or - $n) - 620r (3.12) 


where c = N2/N,. In the previous equations we have absorbed the parameter y in 
the definition of the rescaled time 7 = ty/N2 and introducing additionally è = a/Y, 
B= Bly, 6 = 0/761 = 81/7, th = m/Y da = 82/7 Ne = M2/7- 

Having specified the mathematical context that defines the problem at hand, we can 
proceed following two different, tough complementary, approaches. On the one hand 
one can simulate the full stochastic process, via the celebrated Gillespie algorithm. On 
the other, analytical tools can be invoked to grasp the essential traits of the underlying 
dynamics. 

Stochastic effects fade off for an infinite population amount and the system of equa- 
tions (3.1)-(3.8) tends to its deterministic mean-field solution, as remarked above. For 
a moderate number of individuals, the intrinsic noise due to the discrete nature of the 
individuals can instead play an important role. In some cases, indeed, the stochastic 
process may results in large persistent fluctuations which, with reference to this specific 
case, can drive significant repercussion on e.g. the effect of a chosen pharmacologi- 
cal therapy. As we stated in the introductory section, cycles do emerge at the level of 
the electric brain activity, as registered via the ERPs for patients under pharmacological 
treatment. The claim substantiated in [66] is that such oscillatory behaviors might ap- 
pear due to finite size corrections to the mean-field dynamics, being hence intimately 
related to the intrinsic molecular granularity. Different qualitative patterns might be 
associated to peculiar chemical properties of the drugs, an information which is here 
incorporated into the assigned rates. 

The onset of these oscillations depends on the nature of the equilibrium points of the 
mean-field equations. Finite size induced fluctuations may give rise to persistent cycles 
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when the equilibrium point is a stable focus: Intrinsic noise perturbs the systems, and 
prevents it from approaching its natural asymptotic solution, via damped oscillations. 
Under this condition, a resonance may develop which yields to oscillations at a given, 
dynamically selected, frequency*. Aiming at identifying the necessary condition for 
these oscillatory regimes to develop about the mean-field equations (3.10)-(3.12), one 
needs to single out the parameters’ values which correspond to complex eigenvalues of 
the Jacobian matrix associated to the mean-field systems. In our case the elements of 
the Jacobian matrix reads 


lG) = åC- phr) - Bom, —(1+6)b%-b1- th 

ha) = db, + bA- pr- dî) 

Ji3(@") = Bop, —(1+ 6)b>- th 

Ial) = AAL- prr) + Êk] 

Io(d) = cad + ÊU- pr- en) 

Jo3(@") = BOR, 

Jal) = Ph- te 

Jso(6*) = 0 

Iss(0) = Ph- h- by (3.13) 


where d* = (44, Phr, G4) denotes the equilibrium point of system (3.10)-(3.12). To 
analytically addressing the full problem proves rather cumbersome. For this reason, 
in the next section, we start by discussing a simpler case study, and eventually move 
towards the general settings. 


3.2 On the case with 72 = 0: Neglecting the inward migration of inactive molecules 


Molecules of species H can appear the region of the interaction, as resulting from 
two independent mechanisms: Diffusion may occur as prescribed by equation (3.8) as 
well as chemical inactivation (3.4). Imagine this latter to dominate over the former, an 
assumption which we here exemplifies by setting to zero the migration contribution, i.e. 
n2 = 0. This working ansatz enables us to write down closed analytical expressions for 
the two equilibrium points d and di A straightforward algebraic manipulation leads 
to: 


CA = (tr) 
ô+ 9d, + AM 
A[d2(1 +6) + | auiii 


i ui EEE: x 
b | A[do(1+6)4+ 71] + G[(14+6)(1- 62) +61] (1+6) +h 


4 The importance of this, rather general, property is being addressed in [6, 7] and will be touched upon in 
the forthcoming discussion. 
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As previously mentioned, a stability analysis represents the first mandatory step to 
identifying the parameters range which yields to the oscillations. As concerns di one 
can easily calculate the three eigenvalues A; (2 = 1,2,3) of J (07) which respectively 
reads: 


en 
fji + 01 
ò âñı +854 
\ 7 — [i +5, + Base td |a VDI 
233..° > "Se eee ee 


2 


where 
G28? (+1) , 2001 - c) 
(Ami + Bb)? din + Bb 
Eigenvalues Az and A3 are always real and negative, so that è is globally stable if 
A, < 0, a condition which immediately translates into: 


Da = (ed) |i do], + 2câĝ 


— (3.14) 
nı at ôi 
The stability of the second equilibrium point d, cannot be handled analytically and 


it is hence carried out numerically. To this end we notice that the characteristic polyno- 
mial associated to the Jacobian matrix (3.13) is a cubic equation of the type: 


aà3+bA° +cA\ + d=0 (3.15) 


a, b, cbeing implicit function of the relevant chemical parameters. One can hence intro- 
duce the so called discriminant A defined as A = -46° d+b?c°-4ac3+18abcd-27a?d?. 
It is well known that the cubic equation admits a pair of complex conjugate roots if 
A < 0. One can then scan over the available parameters space, looking those combi- 
nation which implies a negative discriminant, namely the potential signature of an os- 
cillatory stochastic dynamics. In the next subsections we shall report on our numerical 
and analytical findings. 


3.2.1 Considering the case A = fı 


In the proposed model, chemical reactions (3.5) and (3.6) stem for the local diffusion 
from the inspected target region. As a first working hypothesis, one can equate incoming 
and outgoing contributions and hence assume the simplifying setting where 6) = 7. In 
practice this amount to sampling the interesting parameter region along a preferential 
direction. 

Under this assumption, we set down to study the sign of the discriminant A associ- 


ated to di when scanning the plane (di, by) for fixed N,, No, à, B and 6 (values as in 


the caption of Fig. 3.2). The colored region shown in Fig. 3.2 corresponds to pair (51, 62) 
which yields to complex eigenvalues of the Jacobian J( CA ). The colorcode refers to the 


(absolute) values ofthe corresponding complex eigenvalues. A visual inspection allows 
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Figure 3.2: We report a projection of the admissible parameter space: Here d1 and da are allowed to 
vary. 71 also changes under the constraint 71 = 61.Th e other parameters are keptfi xed: â = 0.9, 
B= 0.6, è = 0.5, th = 61 and f2 = 0.Th e number of simulated molecules is Ni = 4000, N2 = 9000. 
Points falling in the colored region, are characterized by complex eigenvalues associated to the stable 
equilibrium point $*. This implies an oscillatory (damped) convergence tofi xed point. Colorcoding 
measures the magnitude of the absolute value of the eigenvalues’ imaginary part.Th e insert shows an 
enlargement of the main panel to fully appreciate the region for ôi > 1. A critical 65""* exists where 
complex eigenvalues are found for any 6, amount.Th e black line represents the theoretical prediction 
(3.17) which accurately locates the singularity in the parameter space. 


us to highlight two interesting features. First, we notice that for ô> > 0.5, no values of 


è, are found that satisfy the sought condition A < 0. This conclusion, reached on the 
basis of a purely numerical evidence, is also confirmed by a simple analytical argument 


which is developed in the following. Consider the condition 6, << 1, and so focus on 
the leftmost region of the plane in Fig.3. 2. Performing a power series expansion of A 


in 6), and retaining only first order terms, one immediately obtains: 


3? (by - 1)? + cà? 62 + a6 bo - 1)d2 i Sg 5 
jo iii aah, cigs 
[96 - 1) - ad] 


An exact calculation allows one to prove that A is negative for ô> < 0.5. Notice that 
this conclusion holds for any chosen parameters’ setting and has to be regarded as a 


general property of the model at hand”. It is worth emphasizing that condition 62 > 0.5 
also arises when looking for positive real eigenvalues associated to di: In other words 
dy = 0.5 sets the frontier for the global stability of @j, as it follows by inserting 7), = by 
in equation (3.14). 


5 Numerical simulations show that this result holds also when fì + dì. More specifically, there are no 
complex eigenvalues for bo >t If + 61). 


48 Francesca Di Patti 


The second aspect which is brought to our attention when analyzing the scenario 
depicted in Fig.3. 2 has to do with the pronounced peak which emerges at a specific, 


critical value of 52. Let us label such a value with dgrit and notice that it corresponds 
to requiring A < 0 Vò,. This latter observation results in a straightforward strategy for 
analytically evaluating oo We start by rearranging the expression of A as a fractional 
polynomial of 61. In formulae we get: 


14 Ci 

Da Anum in Vico adi 
E 10 Sk 
Aden Vk=o brô 


where the coefficients a; and by are defined as functions of the chemical parameters of 
the model®. 


Aiming at calculating ie one can consider ôi large. Under this assumption, Anum ~ 
140}* + a13615, which can be cast in the explicit form: 


Ain 


$13 
ôi 


x -4(a4 a 1+2(1+cé4 8). | [ 20? 6°53 


+2037 da (46. 7 bt â( 24 3cô2) 80 4 6526) 
+å (51 (1 - 252 - 2càd) + ô + 263(-10 - 8cà + Pâ? 


~ 86 - 6c&G) — 26(-5 + cà — 36 4 2câô)) 
+ Ê (2cà?3,(3cò, 2) -( 265 = 1) (4 + ô + 2ô2 + 5ĉ) 
4â(-1 + 2cô2 + d3(2 + 4c + cò, 4 5cô)))| 


As clearly displayed in Fig. 3.2, the boundaries of the two curves which enclose the do- 
main for A > 0, tend to merge asymptotically (seri > oo) and approach the horizontal 
line that goes through derit, As the aforementioned boundaries are found by imposing 
the discriminant to zero, we argue that gere can be estimated by setting A = 0 in the 
limit for large 6, amount. This fact is also confirmed by looking at the inset of Fig. 3.2: 


As ô is increased the imaginary part ofthe eigenvalues progressively reduces (the color 
contrast fades), which again implies A + 0. In formulae, equating expression (3.16) to 


é The explicit expressions of coefficients a; and by has been worked out with Mathematica and involve 
hundreds of terms. It is hence not explicitly reported in the body of the paper but can be supplied upon 
request. 
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zero and solving for by yields: 
(G+ B)(-1 + 20, + 2câð> + 2C32) 
x (46 448 + 1006, - 200259 — 6Bd5 — AB, 
-16cd do — 4câ? Bd, — 1468755 - Acà 8° da — 20863 
-16câ? f2 + 2024362 — 4362 - 8câ b2 + 6c?42 Bd? 
+8072 + 6c 2a 6262 + 2c? 6362 + AG +586 + 64626 
-4câ? 6,6 — 108606 — 20003526 — 1668” 526 
= 160626 - 12câ?$2ô + 12cĝ?62ô ) (3.16) 


Large value of ôi, i.e. the asymptotic condition we eventually wish to recover, are ob- 
tained by setting the above denominator to zero, which immediately implies: 
ceri 1 
e (3.17) 
2 [1 + c(à + d] 


The above prediction is plotted in Fig. 3.2 and proves accurate in interpreting the nu- 
merics. As afi nal remark, we stress that, though calculated with reference to a specific 
selection of the involved parameters, the patterns displayed in Fig.3. 2 are typical for 
the model under scrutiny. In particular, the dependence of equation (3.17) on both â 
and 8 has been tested in a dedicated series of simulations which confirmed its excellent 
predictive adequacy. 


3.2.2 Mean field vs. stochastic numerical computations 


Numerical simulations of both the the mean-field equations and the original sto- 
chastic model are performed and compared to the scenario depicted above. 

Figure3. 3 shows a typical simulation for a set of parameters which is predicted to 
yield complex eigenvalues of J (CA ). Dashed lines in panels (a)-(c) represent the nu- 
merical solution of system (3.10)- (a 12): The mean-field equations exhibit damped os- 
cillations when approaching the equilibrium point o>. >: Panel (d) reports a typical mean- 
field trajectory projected in the (dr, drr) plane: È, is indeed a stable focus, as one can 
appreciate visually. 

Alternatively, one can perform a full stochastic simulation based on chemical equa- 
tions (3.1)-(3.8). This is achieved through the Gillespie algorithm introduced in ap- 
pendix A, which is equivalent to solving the underlying master equation (3.9). At vari- 
ance with the mean-field system, which formally holds for N + oo, stochastic sim- 
ulations (continuous lines in Fig.3.3 ) display large and persistent oscillations about 
the stationary state. These oscillations emerge due to the discreteness of the simulated 
medium and are hence a typical finite size effect, which is deliberately forgotten in the 
realm of deterministic approaches. 

The qualitative comparison between mean-field and full stochastic solutions points 
to the importance of properly accounting for the effect of granularity in any sensible 
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Figure 3.3: Panels (a)-(c) represent a typical solution of the mean-field system (dashed blue line) 
and the stochastic simulation (continuous magenta line). Panel (d) shows the projection of the above 
solution onto the plane ($T, rT). Parameters used in the simulations are N; = 4000, N2 = 9000, 
& = 0.9, 8 = 0.6, = 0.5, 61 = n = 0.01 and ô> = 0.2. 


biomedical model. Traditional pharmacology, in fact, takes the mean-field viewpoint 
and cannot reproduce those collective behavior that might develop as an instability 
driven by the stochastic component to the dynamics. As imagined in our model, os- 
cillations in the number of Rr molecules reflect in a modulation of the experienced 
perception of pain (and of the associated electrical indicators). 


As already stated, the request for a complex eigenvalue of the linearized dynamics 
constitute a necessary condition for the existence of stochastic induced cycle. The con- 
trary is not a priori guaranteed and dedicated computation tools are to be put in place so 
to access an estimate of the subset of parameters which might drive the resonant effect. 
This is for instance the case of the finite N perturbative expansion pioneered by Van 
Kampen. In the following section we revisit the main step of the derivation (already 
discussed in [66]), and exploit the conclusion to gain some insight into the role played 
by the parameters. 
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3.2.3 Parameters for the resonance condition:Th e large N expansion 


The conditions for the resonant instability reported in Fig. 3.3 are investigated via the 
perturbative expansion, termed after Van Kampen [64, 67], which applies to moderate 
value N of the population ensemble. 

The method relies on the introduction of new set of stochastic (and continuous) 
variables, namely £ = (Er, £r,, ÊRp )> related to their discrete homologous via 


No@r(t) + No T 
tr, = N@rr(t)+VMéry 
na = Nobu(t)+VNo€éx 


A new probability distribution II(€,7) = P(n, t) can be hence defined which obeys to 
the master equation (3.9) now re-written in term of the new variables. In doing this we 
bring into the equation an explicit dependence on the amount of molecules (resp. Ny 
and N2). This fact enables us to perform a perturbative analysis inspired to the seminal 
work by Van Kampen [64], where N7! (Eq. N31) plays the role of a small parameter. At 
the lowest order of approximation, one recovers the deterministic mean-field equations 
(3.10), see [66]. 

The first technical point of the van Kampen expansion concerns the introduction of 
the so-called shift operators, E7, ER} , EZ which obey: 


DI 
9 
Il 


Er f(n, t) = f(nrt1,nerr, ny) 
ER, f(n, t) = f(nr, ner +1,nx) 
Et f(n,t) = f(nr, ner, ny +1). 


The master equation (3.9) is hence cast in the form: 
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The advantage of using the shift operators relies in that they admit a simple expansion 
in the limit for N, (resp. N2) large: 


o, Gee a 


UT = 3 2 Ep a Na dE (3.19) 
= o 1 o 

ztl = 1 x N 1/2 = 1 . 

Ri 1 Den, 2 N; ae, (3.20) 

pilo _ Æ -1/2 _d_ 1 SCAP 

los 1aN; IFA a z+ (3.21) 


Plugging (3.19)-(3.21) into (3.18), after some algebraic manipulation, one recovers at the 
leading order the mean-field equations, formally identical to the ones reported above, 
see equations (3.10)-(3.12). The next-to-leading order result in a Fokker Planck equa- 
tion for the probability distribution II(€, t): 


a Lag (AON) + 5 aL Po agoe 4 (3.22) 
where A,(§) =); M;,;§; (hence linear in €;). The elements of the matrix M are 
My, = â(1 Orr) BOR, (1+6)¢y by a 
My = Plg + BU. - OF - on) 
Mis = -Êr ~ (1+ 6) Or ~ th 
Ma = PIAL- di) + Boke] 
Mx = -clà@r+8(1-61-6)] 


Mo3 = H? Bb, 


Msı = by 
M32 = 0 
M33 = obp- do (3.23) 


while the noise correlation matrix B is symmetric and it is given by 
Bu = ÂP- prr) + Bip (1- Or- bi) 
+(1+6)oroy + dbp tM- br - or) 
By = -CPIE - Phr) + BORC- Ot- O77) 


By = -rH 

Bə = c[âpr(1 È dry) + BOR, (1 = oT = Pd )] 

Bog = 0 

Bss = Pry + doby (3.24) 


Here we have set £, = Er, £2 = Rp- £3 = En. In principle, by solving equation (3.22) one 
can quantify the deviation from the ideal mean-field dynamics. As we aim at elucidating 
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the emergence of regular cycles induced by stochasticity, it is convenient to map the 
problem into an equivalent setting of Langevin type, which takes the explicit form: 


Si AO + n(7) (3.25) 


where the term 77; is a Gaussian noise with zero mean and with correlation given by 


(mT) (7 ))= Byêlr-r) 


Here ¿i,j = 1,...,3 with 1 = T, 2 = Rr and 3 = H. 

To bring into evidence any possible oscillatory regime, we perform a Fourier anal- 
ysis and calculate the associated power spectrum. Following [7], we take the Fourier 
transform of (3.25) obtaining 


~iw€;(w) = È M;f;(w) + 7j;(w) 


where we have denoted by the tilde the Fourier transform. Introducing the matrix ® = 
—iwI — M, the previous relation may be written as 


Z 84; (wW) E(w) = fiw) (3.26) 
and the Fourier transform of n;(w) has the following correlation function 
(7i(w)@;(w ))= Bi;(27)d(w+w ) 


Equation (3.26) gives _ 
Ei(w) = Yo Hw) 
J 


and averaging the squared modulus of £; (w) we obtain the analytical expression of the 
power spectrum 


P(w) =(| &(w)P) = Y 2 Di; (0)Ba(Di,) w) (3.27) 


where DI (w) = ®,,(-w). 

Bound receptors might be assumed as quantifying the drug compound action: The 
largest the number of bound receptor, here labeled Fr, the more effective the analgesic 
in reducing the pain. This is of course a simplistic picture which does not weight the 
certainly relevant contribution associated to the subsequent cascade of reactions (signal 
transduction) engendered by a successful binding event. As we shall be discussing in a 
separate contribution, the model can be made more realistic by accounting for such an 
additional effect. 

Motivated by the willing of unraveling regular patterns in the population of bound 
receptors, we calculate the power spectrum defined by (3.27) for the case of interest 
i = 2. After some algebraic manipulations, P(w) can be put in the form of a fractional 
polynomial of the type 


2 4 
Co + Caw? + Caw 
P3(w) = —T_____ 3.28 
2(w) do + daw? + dwt + dgw® 28) 
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Figure 3.4: A subset of point (lighter region) is identified which leads to an oscillatory behavior of the 
acquired stochastic series for Rr. Such a region is delimited on the basis of an analytical perturbative 
treatment, see equation (3.27). 


where c; and d; depends on the parameters of the model and on the equilibrium point 
O 7. This is a symmetric function with respect to the axis w = 0, and goes to zero 
when w — oo. At most two maxima (peaks) are possible, depending on the number 


of real roots admitted by its derivative. Scanning the parameter space of Fig. 3.2 we 
discover that only a subset of pairs (61, 62) yielding to damped oscillations in the mean 
field regime, do admit a peak in P;(w). The (in principle possible) two peaked profile 
is not being found within the explored parameters range.Th e lighter region in Fig. 3.4 
highlights the portion of space corresponding to the predicted oscillatory behavior for 
Rr. 


To validate our finding we turn to direct, Gillespie based, stochastic simulations and 
choose the parameters that are expected to lead to a peaked power spectrum. This latter 
is calculated averaging over many independent stochastic runs and then compared to 
the theoretical curve (3.27), in Fig.3. 5. A good agreement is found. The existence of 
a pick in the power spectrum implies that the concentration of bound receptors will 
undergo significant time oscillations induced by the grainess intrinsic to the medium 
being simulated. 


3.3. On the general case where f2 + 0 
Let us consider the general scenario where renewing of H molecules is also allowed, 


implying 72 + 0. In this case the mean-field system (3.10)-(3.12) admits only one equi- 


7 Also in this case we do not report the explicit expressions of the coefficients c; and d; as they do involve 
hundreds of terms. 
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Figure 3.5: A plot of the power spectrum relative to species Rr as a function of frequency w. The 
symbols represent the profile calculated from 500 independent runs of the stochastic simulations, 
while the solid line stands for the theoretical prediction based on equation (3.27). Parameters are 
those listed in the caption of Fig. 3.3. 
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where 


D, = DIS fd. Ôi Ĥa + 


A24 4ñ2ô1 [a + ô) (êz +2) + in| 


and A = (1 t ô) t hh 6169 1102 Ore. 

Also in this case we are interested in studying the emergence of cycles and so we 
follow the procedure outlined above. Due to technical difficulties in carrying out alge- 
braic manipulations within the generalized formulation, we proceed with a numerically 
assisted investigation. It shall be however emphasized that the van Kampen calculation 
can be straightforwardly extended to the case where 72 + 0. This case, in fact, is imme- 
diately recovered by replacing M3; > M3) — 2 and M33 + M33 — 72 in matrix (3.23), 
and B33 > B33 + fo(1 - o> — př) in matrix 3.24. Analytical prediction for the peak 
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Figure 3.6: Scanning the parameters space (51, 62): Dark patches are associated to complex eigenvalues 
for ¢*, while a lighter color is used to trace the subdomain where a peaked power spectrum is found 
for the Rr species. Here Ni = 4000, N2 = 9000, & = 0.9, Ê = 0.6, & = 0.5, 7f1 = d1, è = 0.000001 
for panel (a), 7/2 = 0.001 for panel (b), and 7j2 = 0.0099 for panel (c). We notice that also in this case 
the numerically calculated power spectrum is correctly interpolated by the theoretical curve obtained 
within the van Kampen expansion (data not shown). 


profile are hence obtained and here used to discriminate the parameters’ values which 
lead to an amplification of the stochastic oscillations. The results of this study are dis- 
cussed in the following: First we start by considering, 7}; = by, in analogy with the above 
analysis. Then we relax this latter condition and allow 7, f2, ôi, bo to change freely. 


3.3.1 Assuming fı = ôi 


To make contact with the analysis of Section 3.2.1 we start by imposing fı = 6,. Us- 
ing the explicit van Kampen prediction for the power spectrum profile, we scanned the 
parameter plan (61, 62), looking for the presence of a peaked profile. Results are dis- 
played in 3.6 where three different values of the 1), are being considered. Dark patches 
identify region where complex eigenvalues are found for the dynamics linearized around 
*. The superimposed lighter domains define instead the parameters’ range where a 
peak is observed for the power spectrum relative to the molecules of type Rr . 

For small values of 7/2, the obtained diagrams resembles pretty closely that reported 
in Fig. 3.6. Increasing the value of t2, the region which yields to complex eigenvalues 
expands, while the subdomain corresponding to a cyclic dynamics shrinks. 


3.3.2 General case 


As anticipated we now consider the general setting where the parameters 7), 72, di; 


ô can be independently adjusted. Results are reported in Fig. 3.7, where now a three 
dimensional plot is introduced. The other reaction parameters are frozen to the val- 
ues employed before.Th e qualitative conclusions reached holds, however, in general, as 
confirmed by a campaign of simulations, not reported here. A significant, densely con- 
nected, region of the parameter space is systematically detected for which oscillations 
induced by the stochastic component of the dynamics are present. 

Panels (d)-(f) represents the projections of the three-dimensional plots on the plane 
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Figure 3.7: Analyzing the three-dimensional parameters space (1 60, tp). Dark domain is associated 
to complex eigenvalues for d*, while the lighter region represents the subset where oscillations in the 
Rr time series are present. Here Nı = 4000, Nə = 9000, â = 0.9, B 0.6, 6 = 0.5, HL = b1, 
ta = 0.000001 for panel (a) and (b), 172 = 0.001 for panel (c) and (d), 7j2 = 0.0099 for panel (e) and 
(f). 


è, î). Oscillations in the number of bound receptors are frequently occurring when 


Ĥı > 61, a condition that corresponds to having the incoming rate of active species, 
larger than the associated removal probability. In other words, translating this conclu- 
sion into a sound medical interpretation, we might be inclined to expect large oscilla- 
tions in the experienced perception of pain when exceedingly massive doses of analgesic 
are administered. To illustrate the phenomenon, we report in Fig. 3.8 a typical stochastic 
simulation, together with the corresponding mean-field solution. 

Indeed, the above discussion applies to relatively large population amount. When 
instead N2 gets smaller (the threshold being sensitive to the actual value of 72), the 
system approaches a so called absorbing boundary, and the van Kampen ansatz breaks 
down. More specifically, the system switches between the two attractive states of the dy- 
namics (the stationary state d°, and the equilibrium point di characteristic ofthe model 
for 7) = 0), the number of H molecule becoming occasionally zero. This process is re- 
sponsible for continuous, though unpredictable swing in the number of bound receptor, 
yet another dynamical regimes where analgesic might prove substantially ineffective. 


58 Francesca Di Patti 


L L L L L L 
0 2000 4000 6000 8000 “0 2000 4000 6000 8000 
Time Time 
0,05 0,8 
0,7} 4 
0,04 4 
0,6} 4 
0,03 4 
= M0 5- 4 
e e 
o2) dal: | 
‘III Lil || LA 
"| ‘IN nM I LUNI) NI TA ost J 
| 02 fi fi 
% 2000 4000 6000 8000 o Ol D 93 04 
Time T 


Figure 3.8: Comparison of stochastic simulations (solid lines) and the mean-field numerical solution 


(dashed line). Parameters used for the simulations are Ni = 1000, N2 = 9000, & = 0.1, B = 1.67-10°?, 
è = 0.2, f = 6.67 - 1073, f2 = 3.33 - 1079,d1 = 3.33 - 1073 and do = 0.27. 


We notice that this is a general property of the model, which can be encountered also 
in the simplified settings discussed above. To report on this intriguing phenomenon we 
represent in Fig. 3.10 two typical realizations (resp. for 17, + 6, and 7}, = 6, ) for which 
it is found to occur. 

As a final remark, let us stress that the process of signal transduction can be further 
detailed by elaborating on the receptors activity, following the lines of [68]. It can be 
shown (data not displayed) that the main conclusions here reached still hold under such 
general, and more realistic, conditions. 
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Figure 3.9: Comparison of stochastic simulations (solid lines) and the mean-field numerical solution 


(dashed line). Parameters used for the simulations are Ni = 2000, N2 
& = 0.2, ra = 6.67 - 1075, f = 3.33 - 1076,6; = 3.33 - 1073 and d2 = 0.27.Th e stochastic simulations 
display a switching between the equilibrium point of the system and the trivial equilibrium point of 


the system for Ĥ2 = 0. In this case, È = (0.667, 0.923, 0). 


2500, â = 0.1, 8 = 1.67-107?, 
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Figure 3.10: Comparison of stochastic simulations (solid lines) and the mean-field numerical solution 
(dashed line). Parameters used for the RC are Nı = = 1000, N2 = 9000, â = 0.1, B= 1.67-10°?, 
ô = 0.2, 71 = 61 = 3.33- 1073, ro = 3.33 - 1077 and 62 = 0.27.Th e stochastic simulations display a 
switching between the equilibrium point of the system and the trivial equilibrium point of the system 
for 72 = 0. In this case, ø} = (0.5, 0.86, 0). 


Chapter 4 
A stochastic approach to the coupled parent drug and metabolite 
dynamics 


As outlined in the previous chapter, stochastic effects might be important when in- 
terested at describing the interaction between molecules. In particular, we have so far 
analyzed the interplay between analgesics and the corresponding target receptors. In 
this chapter we take one loop forward by developing a model which includes the possi- 
bility for the analgesic to undergo a biotransformation into metabolites, a distinct family 
of molecules which competes with the parent drug for the same receptors. More specif- 
ically, it is important to account for the peculiar role of different metabolites, which 
can arise from the genetic variability of the enzymes responsible for the biotransfor- 
mation. These variations may substantially affect the individual response to the ther- 
apy, as commonly experienced in the medical practice. Detecting genetic variations in 
drug-metabolizing enzymes becomes e.g. essential for identifying individuals for which 
adverse drug reactions to standard doses of certain medications are expected. Individ- 
uals carrying cytochrome poor metabolizer variants exhibit different pharmacokinet- 
ics properties as compared to control subjects. As a result, non-conventional doses of 
medications are to be eventually required so to sustain the involved cytochrome activity 
for biotransformation. Conversely, medications that are not processed via cytochrome 
biotransformation, can be preferentially selected for those patients with potentially im- 
paired cytochrome metabolic capacity. Although we here make specific reference to the 
case of tramadol, the model here discussed is rather general and can be hence invoked 
within other contexts where metabolization and ligand-receptors interactions do occur. 


In the next section we introduce the stochastic model in term of the associated chem- 
ical equations. The underlying master equation is also specified. In section4. 2, we re- 
cover the mean-field system which formally applies to the limit of infinite microscopic 
constituents. The fixed points of the mean-field model are studied, as well as their as- 
sociated stability. As we shall pinpoint in the following, depending on the chemical pa- 
rameters the drug act with a different degree of effectiveness, that we here quantify. Also, 
the transient dynamics present intriguing features, that we bring into evidence. Section 
4.3 is devoted to investigating the role of fluctuations which are analytically studied via 
the van Kampen’s expansion already applied in the previous chapter. Numerical simu- 
lations are performed to corroborate our findings. 
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4.1 Description of the model 


The bloodstream in the vicinity of the target receptor is assumed to be populated by 
two species of molecules, the parent drug tramadol and the main metabolite M1. Notice 
that, in general, metabolization and binding events occur in different parts of the body. 
Diffusion between sequentially ordered compartments should be in principle consid- 
ered. In the present formulation space is however not explicitly incorporated and the 
reactions happen, according to their associated probability, within a unique bulk where 
molecules are uniformly stirred. As anticipated the solely biological processes here ad- 
dressed are hence the metabolization and the reversible chemical reactions between the 
molecules and the free target receptors. As a side comment, we also emphasize that 
competition with other molecular entities dispersed in the medium could be possibly 
included in the picture. This important aspects are discussed in [66, 69], as well as the 
preceding chapter. 

The process of biotransformation through the cytochromes gives rise to the metabo- 
lites. Here we hypothesize that the cytochromes are present in great quantity in the body, 
so that metabolization does not depend on their associated concentration, and proceeds 
as a spontaneous transformation at constant rate. Denoting with T the molecule of tra- 
madol, and assuming M to label the metabolite of type MI, the process of metaboliza- 
tion is re-conducted to the following chemical reaction 


T&M 


where a is the reaction constant. This is the parameter which quantifies the ability of 
the body to metabolize the drug and can be hence supposed to be intimately connected 
to the genetic polymorphisms of the cytochrome CYP2D6. 

Tramadol sailing in the bloodstream can eventually encounter a free target receptor, 
hereafter labeled Rp. Following a successful binding event the receptor Rp changes 
into an occupied element Rr. In formulae: 


T+ Rr & Rr 


We further assume that during the binding the chemical properties of T remain un- 
altered. It may hence occasionally occur that Rr undergoes the inverse transformation 
by realizing a, still active, 7 molecule. This assumption translates into 


Rr T+ Rp 


The two parameters (3; and 7; are the constant reaction rates. 

To complete the model, we have to include the interaction between metabolites and 
receptors. The medical literature reports on the specificity of tramadol and metabo- 
lites to the different classes of receptors involved in pain mechanism, and their role in 
achieving analgesia. Here, we set down to consider a simplified scenario, where the par- 
ent drug and metabolite bind to the same type of receptors. More specifically, we shall 
assume, in analogy with the above, the following reactions’ scheme for M: 


M + Rr “Ry 
Ru M+ Rp 
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Figure 4.1: Chemical equations and schematic representation of the model. A tramadol molecule (T) 
turns into metabolite (M) with rate a, and it can bind to a free receptor Rp with rate 61.Th e sponta- 
neous detachment of the compound Rr occurs with rate y1. Molecules of type M bind and unbind 
to Rr with rate 82 and q2. 


where now Rm is the compound receptor-metabolite, and ( and 2 stand for the as- 
sociated (forward and backward) rates. 


The cartoon in Fig. 4.1 depicts the reactions network that we imagine to characterize 
our model. 


How to quantify the sensation of pain within the proposed picture? As already re- 
marked in chapter 3, we imagine that the more the bound receptors the less the ex- 
perienced pain. In other words, the ideal condition where all available receptors were 
screened by a pool of injected drug (T' or M) molecules, would correspond to achieving 
complete analgesia. In future perspective the model could be complemented by accom- 
modating for the signal transduction pathway, and so accurately representing the neural 
activity steps involved in the process. 


Another comment is mandatory at this point. We have in fact deliberately decided 
to disregard the elimination of tramadol and metabolites from the body. Elimination is 
indeed crucial and leads to thefi nal absorbing state where the concentration of T' and 
M are both zero. However, and being at present interested with elucidating the local 
interaction of drugs and receptors, we hypothesized the elimination to proceed on a 
different (sensibly longer) time scale. Under this working assumption we do imagine to 
focus on a sequence of snapshots of the (relatively faster) interaction dynamics, where 
the global number of microscopic actors can be assumed as constant. 


Moreover, we can certainly assume that the total number of receptors does not change 
with time (degradation is also happening with a different characteristic time). Denoting 
with n; the number of molecules belonging to the i-th species, for î = T, M, Rr, Rm, 
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Rpr, the following constrains are in conclusion put forward: 


Np t+ny t+ NR +NnRm = Ni 
NRE + NR +NRm = No 


where NV; represents the total number of molecules, while N3 refers to the total number 
of receptors. We can use these relations to express nz,, and Nr, in terms of the other 
independent variables, namely nr,, = Ni-nr-nu - Nr, and nr, = No- Ni + 
ny + ny, so that the state of the system is given by the three dimensional vector n = 
(nr, NM, Noy). 

Within this framework, we are able to write the transition probabilities for the sys- 
tem to go from initial state n to the final (allowed) state n/. Such a probability is la- 
beled 7(n/|n). In our system only transitions from n to (nr -1,mnm+1l,nr,) (nT +t 
l, nm, Npr, F 1) and (nr, ny + 1,nR,) can take place. The corresponding nonzero 
T (n\n) entries are 


T(nr-1,nm+1,nr,|n) = anr 
No- N 
T(nr-lnmnar + Ila) = Bing AE or + me) 
Nə- N 
T(nr,nm-1,nerh) = Pony MaMa nr + a) 
T(nr+1,nm,ng,- ln) = npr, 
T(nr,nm+1,ner|e) = 2%(Ni-nr-nm-nrr) 


where N = N; + No. 

Transition probabilities allow us to write down a master equation which governs the 
time evolution of the probability P(n, t), namely the probability of having the system 
in state n at time t. The rate of change of P(n, t) is simply given by the sum of the 
transitions towards n, minus the outward transitions propagating from that state. In 
mathematical notation: 


Í Pint) = T(n|nr + l,nm ne l,nrr)P(nr+ l,nm E l,nrpt) 
+T (n|nr + l,nm,Nrp _ 1)P(np + 1, NM, NRr = 1,t) 
+T(njnr,nm + l,nr,)P(nr,nm + 1, nprr,t) 
+T(n|jnr-1,nm, nr, +1)P(nr-1,nm, npr, +1,t) 
+T (nnr, nm E l,nr,)P(nr,nm 7 1,nRr,t) 


[Tr 1,nu t l,nr,|n) + T(nr = l,nm,NRp + 1jn) 
+T(nr,nm-1,nr,|M)+T(nr+1,nm,nr,-1|2) 


+T(nr,nm +1, nein) |PC, t) (4.1) 


with null initial and boundary conditions. 
We have by now formulated our discrete stochastic model and specified the tran- 
sition probabilities between the admissible states. The (exact) master equation could 
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be, in principle, solved to obtain a closed expression for the expected probability at 
time t. This task turns out impossible and one has to resort to approximate solution. 
In first place, as discussed in the next section, the mean-field limit (namely N} + oo, 
N2 + 00) is recovered following the standard procedure applied in the previous chapter. 
Then fluctuations around the mean-field dynamics are considered via the perturbative 
techniques. 


4.2 The deterministic limit 


Multiplying both sides of equation (4.1) by nr and summing over all states, we obtain 


d 
di ) nrP(n,t) = bs [Tr T 1, NM, NRyr ni 1jn) 


T(nr L,nm l,nr, In) = T(nr = 1, nm, Rr + Un) |P(n,t) 


where the summation variables have been shifted to simplify the expression. Substitut- 
ing in for the transition rates and remembering that by definition 


)nrP(n,t)= (nr) 


we have 
d 
dt 


Applying the same method to the two other variables, we obtain the following differen- 
tial equations 


(nr) =- (nr) - Bing Ns Ni +nrt+nm))+ 2y1(nR,) - (4.2) 


€ (nw) =a(nr) - P n (No -Ni +nr+nym)) 
+ 2y2(Ni - nr - nm — NRr) (4.3) 
Tinm) = P np(N; Ni +nr+nm))-2y:(nr}). (4.4) 


In the limit N + oo we can replace (n;n;) =( n;)(n;) for every i, j in equations (4.2)- 
(4.4). In this way (n;)/N becomes the deterministic variable d;, and we can write the 
mean-field system as: 


Lor = -adbr — Bidr(0 + or + om) + 21r 


È du = adr + 272(~ - pr- $m- br,) — P20m(0+ dr + dm) (4.5) 


È dr = frbr(o+ or + bu) -Mór 


where o = (N - N,)/N and y= N;/N. 
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Figure 4.2: Metabolites and receptor densities as function of time. The solid (red) lines represent the 
stochastic simulation while the dashed (black) lines correspond to the numerical solution of the mean- 
field system. Parameters used are a = 0.12, G1 = 0.23, G2 = 0.33, yı = 0.4, y2 = 0.35, Ni = 4500 
and Nə = 2000. 


Figure4. 2 shows the comparison between the stochastic behavior of the model and 
the mean-field one, as calculated by numerical integration of equations (4.5), dashed 
(black) line. The continuous (red) line is a typical stochastic simulation obtained through 
implementing the Gillespie algorithm (see appendix A). For each species, the two pro- 
files overlap well:Th e approximate mean-field theory and the stochastic simulation 
agree, a part from corrections due to thefi niteness of the simulated medium. As T is 
not continuously administered, its (number) concentration (as well that of Rr) decays 
to zero. Conversely, the densities of the other species settle down to a steady-state value. 
This latter value and its stability properties are calculated in the following, where simple 
speculations on the medical relevance of ourfi ndings are also going to be addressed. 


4.2.1 Analysis of the macroscopic equations 


To find the equilibrium point of the macroscopic equations, we set the time deriva- 
tive to zero in system (4.5) and solve, obtaining the point $* = (0, ©},;0) where 


—(272 + 820) + (272 + B20)? + 8 72/2 


ou = 23, 


The stability of this point can be deduced from the Jacobian matrix 


-a- b10 - Bibi 0 2 
Jb") =] @=2%2- Body 2V2- p20 - 282,04 -2V2 
Bio + BPm 0 2% 
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whose eigenvalues are 


Ay = -y (272 + 20)? + 892982 (4.6) 
(271 + a+ Bio + Bipi) + V (271 +a + B10 + bipu) -8710 
A2,3 = —— do,_— oo ———_.. 47) 


2 


These values are clearly negative and real, proving the system has a globally stable equi- 
librium point. 

As explained in section 4.1, we here assess the effectiveness of the pharmacological 
treatment by measuring the number of bound receptors. However, recalling that the 
initial dose of tramadol is completely metabolized at equilibrium, we shall be solely 
interested with the quantity @};,,. To visualize the asymptotic stage of the evolution we 
refer to the plan (©r,, dry) and therein trace the bisectrix (dashed line in Fig. 4.3a). 
Above the diagonal, hy > Pkp the drug works better and the patient experience less 
pain. Such a condition realizes if 


Bo +492 
2( 2 - 272) 


and 03 > 22. This means that the forward binding rate for the metabolite must be (at 
least) a factor two larger than the corresponding dissociation constant. Moreover, the 
initial dose of administered drug has to be chosen so that N, is larger than (at least) 
N: 2 / 2; 

As it is shown in Fig. 4.3a, the equilibrium point is confined on the line Nok, + 
No, 7 Na = 0, which also defines the sub-domain of the plane which can be visited 
during the transient dynamics. Above that line in fact the positiveness of the variables 
is guaranteed. In Fig. 4.3b we project the numerical solutions of the system (4.5) in 
to the plane RØR. The trajectories of the mean-field equations evolve towards the 
attractor. Starting from an arbitrary initial condition characterized by a generic value of 
NR, and nr,,, the system gets apparently trapped into a transient phase which displays 
an almost constant number of bound receptors, sensibly different from that eventually 
achieved at equilibrium. Indeed, as testified by Fig. 4.3b, the number of bound receptors 
initially shrinks and only after, due to the action of newly injected chemicals, starts 
growing to approach the fixed point. This setting could correspond to mimicking the 
condition where a patient is exposed to a treatment which closely follows a preceding 
drug administration. 

Furthermore, the characteristic time of equilibration can be estimated as the (abso- 
lute value of the) inverse of the maximum eigenvalue (4.6)-(4.7). This is an interesting 
indicator as it quantifies the ability of the system to eventually attain the asymptotic 
condition where the largest number of receptors is occupied. The data in Fig. 4.4 show 
that the relaxation time decreases as the metabolization rate increases, thus suggesting 
that the administered drug acts more rapidly for extensive. 


1 2 


4.3 The van Kampen expansion 


Corrections to the mean-field dynamics can be calculated by resorting, with analogy 
to chapter 3, to the van Kampen’s expansion. The main idea is to write the variables n; 
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Figure 4.3: Panel (a) reports the projection of the equilibrium point into the plane ¢z,,¢R,,-.Th e two 
domains separated by the bisectrix are respectively labeled less/more pain, according to the prescrip- 
tions of the model (see main text) The colored (yellow) region denotes the portion of the plane where 
trajectories are not allowed. The solid black line correspond to the condition N¢z,, + NOR, -N2 = 0 
(here N; = 2000 and N2 = 5000) Th e equilibrium point (black square) belongs to this line. Panel (b) 
represents the projection of the trajectories on the plane Ø Rp r y for different initial conditions. Pa- 
rameters used for the numerical integration of the mean-field system are a = 0.7, (1 = 0.2, G2 = 0.9, 
yı = 0.8, ya = 0.05, Ni = 2500 and No = 2000. Inset:Th e time evolution of the total number 
of bound receptors, , is reported.Th e curve is traced with reference to one specific initial condition, 
namely the red trajectory of the main panel. 


as a sum of two contributions, namely n; = N @;(t) + VN&,, where i = T, M, Rr. 
Here ¢; stands for the deterministic component, while €; relates to the fluctuations. The 
new probability distribution function IT is hence defined as II(£,t) = P(n,t) where 


£= (Er, Em, Err). Moreover: 


To simplify the notation, it is practice to rewrite the master equation (4.1) in terms of 
step operators 


Ee f(n, t) = f(nrt1,nu, ner) 
Bart) = f(nr,nm = l,nrr) 


ER, f (n,t) = f (a7, nm, nRT + 1) A 
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Figure 4.4: The color code (in logarithmic scale) refers to the slowest characteristic time of convergence 
to equilibrium as calculated from the linearized dynamics. We here scan the parameter plan (a, 2). 
The other parameters are chose as 61 = 0.3, y1 = 0.5, y2 = 0.1, Ni = 3000, N2 = 2500 


The master equation can be therefore cast in the form: 


TPt) = (BFE; - NanpP(n,t) 


+( oa Ep. z du - N, +np+nm)P(n,t) 


+E - 12 yr (Na - Na +nr +nm)P(1,t) 
+(E7 E} - crt t) 
5 Ay WON He Be POO (4.8) 


~ 
I 


ptt 


The operators E¥- change n; in n= 1 and so &; in é; +1. They hence admit the following 
representation i in terms of differential operators: 


2 
gži = 1 x= NU? si N° 3 arte (4.9) 


Substituting relation (4.9) into (4.8) and collecting contributions of order VN, one 
recovers the mean-field system of coupled differential equations (4.5). Working out 
the next-to-leading order, namely N, one eventually obtains a Fokker Planck equation 
(FPE) which characterizes the fluctuations around the asymptotic mean-field solution. 
The FPE reads: 

oll o'I 


le) 1 
ao 2 DE (4: (£)II) + 5 D Bo TEIE (4.10) 


where 


AGE) = DM 
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The entries of matrix M, hereafter m;;, those of B, labeled b;;, are function of the 
chemical parameters of the model and read 


mi = -[a+26:@0r +00 + biom] 
mi = -bır 

m3 = -271 

Mə = a- Bodum - 272 

mor = 2020m+020+ Babr + 272 


M23 = -292 

mai = 2bir + bio + biom 

M32 = Bır 

Ms3 = -27 (4.11) 


bi = Bopr+ Bie? + 010rdm + abr + 2% OR, 

b2 = -adr 

big = -[Biodr+Bd1+Bdrdm +2 er, | 

boo = Bop + b20ġbm + Bodmpr + 2729 + (a - 242) br - 2726 ~ 2VPRr 
bos = 0 

bss = Bioro + bibr + birom +2ndrr (4.12) 


Equation (4.10) can be solved explicitly: The obtained probability distribution I(€, t) 
is a Gaussian and it is hence completely specified by itsfi rst and second moments. In the 
next section we shall calculate the associated moments explicitly and test the adequacy 
of the predictions versus direct simulations. 


4.31 Analysis of the fluctuations 


To characterize the moments of the distribution we proceed as follows. We multiply 
both sides of the FPE (4.10) by é; (resp. €;€;) and integrate over all £. One then recovers 


the equations for the mean value of the fluctuations (£;), as well as for the associated 
correlations, (£;É;). 


The evolution of thefi rst moments is found to be governed by the following equa- 
tions: 


Llen) = mir) + mi2(Em) T Mazs(ÊRr) 


Llen) = Malr) + Ma2(EÉm) + mas(Érr) 


Llen) = ms3i(Er) +M32(Em) + mg3(Érp) (4.13) 
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while the the second moments obey to 


Lig) = 2mu(£7) + 2maa(ErEm) + 2713 (Errr) +bu 


T terén) = mai(Ep) + (mar + mo0)(ErÉm)+ Mo3(ErEr,) + M2(Exr) 
+mi3(EmErr) + b12 

Serger) = mgi(67) + Msz(Érém) + (mu + M33 )(ErEre) 
+my2(EuErr) + mis(Ér) + dis 


d 
aa (€) = 2ma(Erêm) + 2m22(Eq,) + 2mM23(EmEr,) + b22 


<éwEer) = msi(Er&u) + mo (ETER) + m32(£%,) 
+(M22 + m33)(EmÉnr) + mas(Ehr) 
£ (Ee) = 2mail(ErEr,) + 2mMz2(EmErr) + 2mM33(ERr) + b33 (4.14) 


where the elements m;j and b;; are the entries of matrices (4.11) and (4.12). 

The above system cannot be solved analytically (indeed we cannot even cast the 
mean-field solution in a closed analytic form). However, being interested in the fluc- 
tuations around the equilibrium point, once the initial transient has damped out, one 
sets to zero the time derivatives in system (4.13) and evaluates the coefficients m;; at 
the equilibrium point ¢”. It turns out that (€7)*’ = (£m) = (Err). Proceeding in 
a similar fashion with system (4.14), one readily finds that all the second moments are 
zero but (£7,)5 which instead reads 


(2) = Bobhy + (B20 — 272) di + 272 
M = 2(2029%, + Boo + 272) 


In this latter case the stationary probability distribution IT()* is given by 


st _ ee È EM 
HO” = Jamey eap] | Fa 


Figure 4.5a shows the projection of the stationary probability distribution II" on 
the plane (rF, ru). As it can be appreciated by visual inspection, the dispersion 
occurs along the direction given by N@&,, + NOR, — N2 = 0 which also contains the 
equilibrium point. In Fig. 4.5b the stationary probability distribution II** is plotted as 
a function of Em. The figure testifies on the predictive ability of equation (4.15) here 
depicted with a solid line, which is shown to interpolate correctly the numerical data 
(symbols). 

Imagine now to partition the plane (rr, Oras) into two regions separated by the 
bisectrix. Moving above the diagonal, the number of screened receptors increases which 
in turn implies reducing the pain, within our simplified scenario. Fluctuations can fa- 
cilitate the road towards analgesia, as outlined in Fig. 4.5a. The stationary probability 
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Figure 4.5: Panel (a) shows the projection of the stationary probability distribution II° on the orig- 
inal plane frrprm-Ih e solid black line represent the bisectrix of the plane. Panel (b) reports the 
stationary probability distribution as function of Éy.Th e solid line represents the theoretical predic- 
tion based on equation (4.15). The symbols refer to the Gillespie like numerical simulation. Here, the 
parameters are a = 0.5, 81 = 0.2, B2 = 0.3, y1 = 0.5, y2 = 0.1, Ni = 3000, and N2 = 1000 


distribution can be be hypothetically employed to quantify the probability of entering 
the region in (rF, rm) where the drug effect is supposedly more pronounced. This 

latter probability corresponds to the area of the distribution above the bisectrix and is 
quantified in p = 0.196 for the chosen parameters’ setting. Interestingly, although the 
mean-field solution predicts a stationary condition characterized by a pronounced sen- 
sation of pain (fp = 0.129 and @hy = 0.121), there is a nonzero probability that, due 
to fluctuations, the system enters a region where pain is partially hindered. Obviously, 
this speculation applies as long as N is finite (though large). This is for instance the case 
where a local stimulus is applies which interests a finite patch of neurons (see for in- 
stance the whisker stimulation experiment [28]). In the general case where the number 
of involved receptors is exceedingly large, the mean-field dynamics takes over and the 
aforementioned distribution shrinks to a delta. 

In the next chapter we shall move to considering a different case study, where the 
stochastic modeling prove fundamental. In particular, as anticipated in the introduc- 
tory section, we will concentrate on analyzing an extended set of autocatalytic reactions 
which can be hypothesized to occur in a minimal model of cell. 


Chapter 5 
Extended auto-catalytic networks 


Autocatalytic reactions have long fascinated physicists and chemists because of their 
unique features [70]. A chemical reaction is called autocatalytic if one of the reaction 
products is itself a catalyst for the chemical reaction. Part of the reason for the interest 
in these types of reactions stems from the fact that even if only a small amount of the 
catalyst is present, the reaction may start off slowly, but will quickly speed up once more 
catalyst is produced. If the reactant is not replaced, the process will again slow down 
producing the typical sigmoid shape for the concentration of the product. All this is for 
a single chemical reaction, but of greater interest is the case of many chemical reactions, 
where one or more reactions produce a catalyst for some of the other reactions. Then 
the whole collection of constituents is called an autocatalytic set [71]. In addition to 
the interesting properties of autocatalytic sets, there is also an intriguing possibility that 
“bootstrap” reactions such as this may have had an important role in producing complex 
or self-replicating molecules required for the origin of life on Earth [72, 73, 74, 75]. 

Theoretical studies of the properties of autocatalytic reactions are typically of two 
kinds. In the first, rate equations for the reactions are written down and these are ei- 
ther solved numerically or their properties investigated using the techniques used in the 
study of dynamical systems. An alternative is to carry out computer simulations of the 
actual reactions themselves. However, as outlined in the preceding discussion, a third 
strategy can be also employed: Using methods from the theory of stochastic processes 
an analytic approach to the full model (and not just the mean field version) is possible. 

Occasionally, oscillatory behaviors manifest for a moderate number of constituents, 
a phenomenon which arise from the discrete nature of the system. As we pointed out 
earlier, these oscillations are purely stochastic in origin. The main tool that is here used 
to elaborate on their emergence in the autocatalytic setting is again the system-size ex- 
pansion of van Kampen [64, 76] (see chapter 3). 

More specifically we shall consider the autocatalytic reaction scheme introduced by 
Togashi and Kaneko [77, 78]. In most autocatalytic reactions there are two types of 
constituent: the autocatalytic and the substrate. The number of the latter type are kept 
constant by continually feeding them in, however the former are not injected nor ex- 
tracted from the system. In this sense the system is closed as far as the autocatalytic 
constituents are concerned, but open for the substrate. In the scheme that Togashi and 
Kaneko investigate, the reactions are cyclic, with k constituents X1,..., Xx reacting 
according to X; + Xi4, > 2Xj41 with X,4, = Xi, i = 1,...,k. The chemicals are 
assumed to be in a container which is well-stirred, but with the possibility of diffusing 
across the surface of the container into a particle reservoir. 
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Togashi and Kaneko [77, 78] used only computer simulation to study their proposed 
reaction scheme. At variance, and following the recipe of chapter 3 we shall proceed 
analytically by writing down the master equation for model at hand, and then studying 
it through a systematic expansion in N~'/?, where N is the system size. To leading 
order one finds the rate equations which appear in [77], and to next-to-leading order 
a Langevin equation which describes the fluctuations about the stable fixed point of 
the rate equations. From previous work we expect that (i) this first-order correction 
will be sufficient to yield results which are in good agreement with simulation data, (ii) 
the large amplitude of the oscillations can be understood as a resonant effect. One of 
the strengths of the technique is that the next-to-leading order corrections give linear 
Langevin equations which can be analyzed exactly for arbitrary values of k. 

As an additional step, we shall extend Togashi and Kaneko model to explicitly in- 
corporate the notion of space. In doing so, the diffusion of the molecular species par- 
ticipating to the dynamics is consistently accounted for and spatio-temporal collective 
motions unraveled via a straightforward application of the van Kampen analysis. Sec- 
tion 5.2 is devoted to discuss this latter topic, combining the analytical and numerical 
viewpoints. 


5.1 Enhanced stochastic oscillations in autocatalytic reactions 
The autocatalytic reaction scheme described previously can be formulated as 


Xi+ Xii “ 2X41 Xea = X1 


X, E i=1,...,k. (5.1) 


Here r;, a; and 0; (with rx+1 = r1), are the rates at which the reactions take place and 
E is the null constituent. Such constituents have to be included so that the number of 
molecules of type X,, n;, are all independent. Figure 5.1 shows a schematic representa- 
tion of these reactions. If the size of the system is denoted by N, then Df, n;j+ng = N, 
where n p is the number of null constituents. While N is fixed, ng may vary as the to- 
tal number of molecules changes with time. In practice, ng does not explicitly appear 
in the formalism; it is always replaced by N - YÈ_, n;. The rate constants a; and 3; in 
equation (5.1) represent the interactions of the system with the particle reservoir outside 
the container. In effect a; and (3; are the rate at which molecules appear and disappear 
from the system, and thus are analogous to birth and death rates. 

As an aside, we note that reaction rates which result from a binary encounter should 
be scaled by the volume of the system, V. That is, r; + r;/V. This follows from a 
straightforward kinetic theory argument [79]. This is an innocent modification as far 
as this study is concerned, since it is carried out at constant volume, but it becomes 
crucially important when the volume is allowed to change, as it does in the analysis of 
the phase transition reported in [77, 78]. 

The state of the system is labeled by the set of integers {n1,...,nx} and, under 
the assumption that the transitions from this state to any other only depends on these 
integers, the system is Markov and may be described in terms of a master equation. In 
constructing the master equation we need to give the transition rates T (n'|n) from the 
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Figure 5.1: Schematic representation of reactions (5.1). 


state n to the to the state n’, where n = (11,..., Nng). If the system is well-stirred, so 
that the probability of a reaction taking place is proportional to its rate and the number 
of reactant molecules, then from equation (5.1) these transition rates are 


Tm... 1, min +1... rale) = ria N ; 


k 
Pinete) (1 = EL). 


T(ni,...,n-1,...,ngm) = bi (5.2) 


N . 
The master equation for the probability that the system is in state n at time t, P(n, t), 
may now be written down: 


k 
ee) = 2 (E aa- 1) [T(n1,...,M;-1,n41+1,...,ne|n)P(n,t)] 
k 
+Y ( 1-1) [T(n1,...,1+1,...,nx|m)P(n,t)] 
i=1 
k 
+D (E 1) [T(11;...,t-1,..., lm) P(n,t)] (5.3) 
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where £#! are again the step-operators as introduced in Chapter 3. Their definition is 
hereafter recalled: 


EE f(n) = f(n,,...,14 + 1,...;Nk). (5.4) 
As previously remarked, equations such as (5.3) are difficult to analyze, but if one is 
particularly interested in large or moderately sized values of N, then the system-size 


expansion provides an elegant way of encapsulating the essential aspects of the model. 
The key assumption of the method is to write [64] 


&(t) 
VN 


= Pi(t) + (5.5) 


From this relation, limy..(n;/N) = @;(t), the fraction of the molecules which are 
of type X; at time t, within the mean-field (N + oo) limit. The fluctuations about 
these are assumed to be Gaussian, hence the 1/V/N in equation (5.5). This assumption 
applies as long as the system evolves reasonably far from the (absorbing) boundaries, so 
that the probability density functions of the X; is Gaussian. In other words, stochastic 
extinctions cannot be captured by our perturbative calculation. 

Substituting equation (5.5) into equation (5.3) allows us to expand the master equa- 
tion as a power series in 1/\/N. We here recall that the step operator can be approxi- 
mated as: 


ii LO, g VO 
E “la ae tan vat (5.6) 


If we set P(n,t) equal to II(É, t), the left-hand side of the master equation becomes 
[64] 


dP(n,t) _ OU(E,t) _ k OU(E,t) dd; 


d ôt VND ae Og; dt Oo 


Applying the ansatz (5.5) to the right-hand side of equation (5.3), the step—operators 
(5.4) take the form (5.6), the n; in the transition rates (5.2) are replaced by @; and €; 
using equation (5.5). This yields the following terms 


(a) Terms of order N!/2; 


le) le) k le) le) 
Z (È = -Qi (1 = xe) OE; + sof} me 
(5.8) 


iI Me 
pei 


(b) Terms of order N~! and involving first order derivatives 
k 


X {ria x E 41 + Tibnin 


i=1 


Paid IEn e + Qi A ad jt Bize zs} me t) (5.9) 


lo, 
Ei — Ti4l Pi 


hi El 


a 
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(c) Terms of order N~! and involving second order derivatives 


1 & 82 82 02 
z 2 [ribo E + TA — 2 | 
sa(1-) 2 o) E de E Bidi sa} Me t). (5.10) 


If we equate terms of the same order in 1/V N on the left- and right-hand sides, 
the leading order gives 


di 


k 
= (Fidi1- Tis Pin) di + Qi ( = 2 è) - Bibi, (5.11) 


where T = t/N is a re-scaled time. This equation is a deterministic equation for the 
fraction of molecules which are of type i. It agrees with that of Togashi and Kaneko 
[77], if one takes into account that their equations are for concentrations and so contain 
the (constant) concentrations of the species in the reservoir. There is also an additional 
term J}; ; in equation (5.11), which is typically present when mean-field equations 
are derived in systems with a fixed size, but not in the phenomenologically postulated 
form. For small concentrations it will not be important, but clearly it will have an effect 
as the ceiling on particle numbers is felt, reducing the number of molecules entering the 
container from the reservoir, as it should. 

The terms of order N~! in equations (5.9) and (5.10), are now identified with the 
remaining term on the right-hand side of equation (5.7). This resulting equation is a 
Fokker-Planck equation: 


oll maul 
Or a Dia DI [Ai (€) II] + Di By EOE; (5.12) 


From equation (5.9) we see that the A;(€) are linear functions of the £; and from equa- 
tion (5.10) that the B;; are independent of them. Explicitly: 


k 
Ai(€) = (rider = riria) Ei + TiQifi-1 — Tei Qiii — & a C= Pic. (5.13) 
and 
—TiGi-1Pi ifj=i-1 


n Tis 1 PiPisr + TiPiPi1 
Pu =] +0; (1- Ths by) + Bibs ifj=i ae 
-Ti+1 iĝi ifj=i+1 


In equations (5.13) and (5.14), @k+1 = @1 and k+1 = 1, which follows from the cyclic 
nature of the model. 
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Since the A;(&) are linear functions of the £; we may write them as 
k 
j=l 


This means that the probability distribution at next-to-leading order, II(€,7), is com- 
pletely determined by the two kx k matrices M and B, whose elements are independent 
of the €;, and only functions of the @;. For our purposes, in complete analogy with the 
preceding discussion, we need to Fourier analyze the fluctuations. It is hence more con- 
venient to characterize the fluctuation in terms of the equivalent Langevin equations: 


k 
Bs 2 Mij&;(T)+ (T), (5.16) 


dr 


where M is a k x k matrix which can be found from equations (5.13) and (5.15), and n; 
is a Gaussian white noise with zero mean and correlator 


(ni(7)n;(7))= Bid (T-T), (5.17) 


and B;, is another k x k matrix given by equation (5.14). 

Equation (5.16), is a stochastic differential equation for the deviation from the mean- 
field behavior. It is the analysis of this equation, together with the mean-field system 
(5.11), that allow us to describe the stochastic aspects of the autocatalytic reactions in a 
quantitative way. 


5.11 On the fluctuations 


In their numerical studies, Togashi and Kaneko [77, 78] looked at the simplest case 
of the model where the rates r;, œ; and 8; were the same for all chemical species. To 
illustrate our method we will do the same, and so from now on we will drop the subscript 
i on these constants, but it should be clear that our analysis also applies to the general 
situation where they are different for each species. With this choice, the deterministic 
equations (5.11) have a single fixed point: 


Pd 


o č a 
_ B+ka' 


where the asterisk denotes the fixed point value. 

In principle the matrices M and B are time dependent, since @; is. However, in 
practice we are just interested in the fluctuations about the stationary state, once the 
initial transient has died out. This is equivalent to replacing ©; with its corresponding 
equilibrium solution. In our case the asymptotic value of @; is specified by equation 
(5.18). Hence, M and B are given by 


(5.18) 


Mo Mı Ma Ma. ma Mg 
m3 Mo Mi Ma M2 Ma 
Ma Mz Mo Mi m2 Ma 
M = m2 M2 Mz Mo... m2 Mal, (5.19) 
m2 M2 M2 Ma Mo Mi 
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where 

Mo = -a- 

m = -a-r¢* 

ma = -Q 

mg = -a+rd¢* (5.20) 
and 

bo di 0 0 bi 

bi bo by 0 0 

B- 0 b bo 0 0 , (5.21) 

0 0 0.. bo bı 

b 0 0.. bi bo 
where 


bo = 2r(d*)° +8 + a(1 - ko") 
bi = ee: (5.22) 


If N is so large that the fluctuations are completely negligible, then the system tends 
towards a state where the fractions of the chemical species in the system are equal, and 
given by equation (5.18), and stays there. Of course, if N is finite this is no longer the 
case and there are fluctuations about this stationary state, and as we will see these can be 
significant even if N is quite large. Since these fluctuations are expected to be oscillatory, 
we begin their analysis by taking the Fourier transform of equation (5.16) to find 


k ~ 
x (-iwd;; - Mij) Elw) = iw), (5.23) 


where the f denotes the Fourier transform of the function f. Defining the matrix 
—iwd;; - Mj; to be ®;;(w), the solution to equation (5.23) is 


El) = E oo). (524) 


To identify the dominant frequency of the oscillatory behavior, we compute the 
power spectrum for the i-th species, P;(w), from equation (5.24): 


> > D7 (w) Ba (8); (7), (5.25) 


j= 


Pw) = (IEP) 


= 

~ 
Il 

a 


Since ® = —iw/ — M, where J is the k x k unit matrix, and since M and B are indepen- 
dent of w, the structure of P;(w) is that of a polynomial of order 2k divided by another 
polynomial of degree 2k. The explicit form of the denominator is | det ® (w)|?. 

From previous investigations of fluctuations of a similar kind [6, 8, 7], we expect that 
the fluctuations about the stationary state (5.18) will be enhanced by a resonant effect: 
For values of w for which | det ® (w)| is a minimum, the power spectra will show peaks 
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which correspond to larger than expected fluctuations at that frequency. This effect was 
first conjectured by Bartlett [80] in the context of the modeling of measles epidemics, 
and later elaborated upon by Nisbet and Gurney [81], who called these stochastically 
induced cycles, quasi-cycles. However, as already emphasized, it is only in the last few 
years that explicit calculations within the system-size expansion have been carried out 
and a quantitative understanding of the phenomenon has emerged [6]. 

To understand the analytic structure of the power spectra, we begin by supposing 
that we can neglect the effects of the numerator on the right-hand side of equation 
(5.25), and simply determine the dominant frequency by looking for the value which 
minimizes | det ® (w)|. The effect of the numerator will be to shift this frequency; we 
are assuming as a first approximation that this shift will be small, as indeed it has been 
found to be in some cases [6]. If \; are the eigenvalues of M, then the denominator of 
the expression for the power spectra may be written as 


| det ® (w)|? = TI (-iw-A;) (iw - AZ) . (5.26) 
j=l 


Since M is real, the A; will be real or come in complex conjugate pairs, so that the 
products in equation (5.26) has one of two forms: 


(i) If A is real, the two factors involving this eigenvalue give (w? + 7). 


(ii) If A is complex: A = AR + 7,7, the four terms involving A and \* give 
laa? 0 i (5.27) 


The resonant effect has its origin in the structure of the factor coming from the com- 
plex eigenvalues shown in the expression (5.27). It is smallest, and so gives the largest 
contribution when it is in the denominator, for frequencies which satisfy 


w =A? - AR- (5.28) 


If there are several pairs of complex eigenvalues and their conjugates, the largest contri- 
bution should come from the pair for which ArAr is smallest. Clearly this will only be 
approximately true since, not only are we neglecting the numerator, but also the factors 
(w? + A?) coming from real eigenvalues, as well as those coming from other complex 
conjugate pairs. However, as we will now see by looking at two specific cases, k = 4 and 
k = 8, these approximations appear to be remarkably good. 

We study the cases k = 4 and k = 8 because they are the smallest even values of 
k for which one complex conjugate pair and two distinct complex conjugate pairs, re- 
spectively, exist (there are two complex conjugate pairs for k = 6, but they are equal, and 
three for k = 8, but two of these are equal). We therefore expect to see one peak in the 
power spectra when k = 4 and two when k = 8. Our analysis, and the accuracy of our 
approximations, can be directly checked by numerical simulation of the chemical reac- 
tion system (5.1) by use of the Gillespie algorithm [79, 82]. This produces realizations of 
the stochastic dynamics which are equivalent to those found from the master equation 
(5.3). Averaging over many of these realizations gives us power spectra after Fourier 
transformation, which are exact to a given numerical accuracy. We now investigate the 
two cases k = 4 and k = 8 in more detail. 
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Figure 5.2: Time evolution of selected species, i = 1,2,3,4 in clockwise order for the case k = 4. 
Here r = 10, a = 8 = 1/64, N = 8192.Th e dashed line indicates the meanfi eld solution. The 
species display a clear oscillatory trend about their meanfi eld values. A paired synchronization, (1,3) 
vs. (2,4) rich states, is also visible, as already observed in [77, 78]. 


Power spectra when k = 4 


The time evolution of the species is depicted in Fig. 5.2. This clearly displays large os- 
cillations which we aim to investigate analytically. Before beginning this analysis, we 
observe that species 1, 3 (odd) and 2, 4 (even) are paired together and move up and 
down from the reference mean-field level in a synchronized fashion.Th is fact was al- 
ready recognized in [77, 78] and shown to drive successive switches between the 1-3 or 
2-4 rich states, close to the absorbing boundary, i.e. when a small number of molecules 
is simulated. The rate at which the changes occur is controlled by the diffusion param- 
eter. However, the details of the transitions stem from a purely dynamical effect which 
cannot be captured within the perturbative analysis developed here. 


Let us now turn to analytically characterizing the aforementioned oscillatory regime. 
To this end we begin by determining the eigenvalues of the M matrix. 


We note that is a circulant matrix [83], and therefore its eigenvalues are given by 


k ş È 
M= ge OO BED (5.29) 


j=l 


where m1; is the element of M in the first row and j-th column. In fact, M is not the 
most general form of circulant matrix; (k — 3) entries in each row are equal (to mz). 
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This leads to a simplified form for the eigenvalues: 


k-2 
i = made +m e" tm > elridtik 
j=2 
ti -2ril]k sin(370/k 
= mo + mı ec? + m e? -m eat 
(5.30) 
where in the last line Z + 0. Putting in the values from equation (5.20) gives 
af -6-ka, if0=0 
dr { -8-2ird* sin(2r0/k), ifl+0. (5.31) 
For k = 4, the eigenvalues are 
do = -8 - da 
At = -8 = 2ird* 
à = -8 
às = Al. (5.32) 


Within the approximations we have discussed, we would expect that there should be a 
single peak in the power spectrum for any one of the chemical species at a frequency 
given by (see equation (5.28)) 
wi = 4r? (g*)? - pole Ao ae 5.33 

c (0) -8 (81 4a)? b (5.33) 
In Fig. 5.3 we show the power spectrum (for the chemical species è = 2) found by 
averaging over 500 realizations from the Gillespie algorithm, together with that found 
from equation (5.25). The good agreement between the simulation results and those 
found from applying the system-size expansion, shows that the method works well for 
N = 5000. The parameters used in this case were r = 10 and a = 8 = 1/64, which 
gives a value of wo % 4 from equation (5.33). From Fig. 5.3 we see this is a surprising 
good estimate for the position of the peak, given the significant frequency dependence 
which we have neglected to obtain the estimate (5.28). 

Another check of the accuracy of these approximations, and so of equation (5.28), 
is to imagine increasing the parameter 8 at fixed r and a, and asking when w? will 
become zero, and so at what frequency will the peak in the power spectra disappear. 
From equation (5.33) we estimate this to be 


b~ = or 8 -vV2ra, (5.34) 
which equals 0.56 for the values of r and a used in Fig. 5.3. Once again this agrees well 
with the full spectrum which predicts the peak to disappear at about the same value. 
As a final check, we measure the position of the peak from a set of simulations run at 
different values of r. Direct measurements (symbols) are compared to the theory (solid 
line) in Fig. 5.4 and are in good quantitative agreement. Again, we recall that adjusting 
the rate r can be equivalently seen as modifying the volume of the system, which is the 
setting investigated in [77, 78]. 
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Figure 5.3: Power spectrum of species i = 2 when k = 4.Th e analytical curve is shown as a solid line 
and the simulation (average over 500 independent realizations) as symbols. Here r = 10, a = 6 = 
10/64, N = 5000. 


Power spectra when k = 8 


From equation (5.31), the eigenvalues of the M matrix are 


do = -B-8a,r 4=-B, 

` = -B-V2ird*, A7= Xi, 

Ag = -B-2ird*, Ae= A}, 

às = -B-V2ird*, As = Aj. (5.35) 


Since there are two distinct complex conjugate pairs we would expect to find two peaks 
in the power spectra, one at w? = 2r?($*)? — B° and the other at w? = 4r?($*)? — 82. 
For small 3, one peak will be at a frequency V2 times the other. We would also expect 
that the peak at lower frequency would be larger than the one at higher frequency, since 
ARA, is smaller for the former. That is, the pole in the power spectra in the complex 
frequency squared plane is nearer to the real axis for the peak at lower frequency, and so 
should have a bigger effect. So, in summary, our approximations indicate that the peaks 
in the power spectra should be given by 


fe GOO 
(G+ 8a) 
4 DD, 
wi, = ES (5.36) 
+80 


with the peak at w = we larger than the one at w = we2. The results of plotting the full 
spectrum found from equation (5.25) and simulation results are shown in Fig. 5.5 for 
r = 200, a = 1.9 and 8 = 2. This corresponds to peaks at w = 31.2 and w = 44.14, 
according to equations (5.36), which once again agrees very well the the results displayed 
in the figure, as does the prediction that the peak nearest the origin should be the largest. 
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Figure 5.4:Th e position of the peak in the power spectrum (for species i = 2 when k = 4) plotted as 
function of the rate constant r. Symbols refers to the stochastic simulations, while the solid line shows 
the analytical prediction. Here a = 8 = 1/64, N = 5000 


To the best of our knowledge, this is thefi rst time that a double-peaked power spectrum 
has been predicted to emerge as a resonant effect, within a van Kampen type of analysis. 


5.2 Ona spatial model of autocatalytic reactions 


We here present a spatial version ofthe autocatalytic model discussed in the previous 
section.Th is work is developed in collaboration with Pietro de Anna and constitute the 
core of his master thesis essay [84] (see also [85]). Our analysis is closely inspired to that 
of Lugo and McKane [86] 


The idea is to introduce a spatial coarse graining at the level of small micro-cells 
(which total to Q) which are supposed to uniformly cover the volume occupied by the 
vescicle. In each microscopic cell autocatalytic reactions as specified by (5.1) (and re- 
called below) do occur. Migration between neighbours cells is allowed, an ingredient 
which in turn amounts to explicitly account for space. As a simplifying ansatz, we im- 
age a periodic geometry and focus on a chain of microscopic cells situated at the fron- 
tier with the external membrane. In the case of al D setting micro-cell 1 is adjacent 
to micro-cell Q. This simplifying hypothesis is put forward so to restore the transla- 
tional invariance, which shall be invoked in the Fourier based treatment developed in 
the end. Physically one might imagine to simulate the transport of material along a tiny 
shell which is positioned close to the outer edge of the vesicle (a sort of spatio-temporal 
zonal flow).Th e shell is then assumed to communicate with the outside and the inside 
via an effective diffusion mechanism which, however, does not embed space explicitly. 
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Figure 5.5: Power spectrum of the time series for species i = 2 when k = 8.Th e analytical result (solid 
line) is superimposed onto the simulations (symbols), averaged over 500 independent realizations. 
Here r = 200, a = 1.9, 8 = 2, N = 7000. 


5.2.1 The perturbative expansion 


We start by consider, in complete a analogy with the aforementioned Togashi and 
Kanenko model, k chemical species. The autocatalytic reactions for the X} species are: 
Xi + Xi, > 2Xi 
where X{,, = Xi. The index i labels the cell where the reaction is supposed to occur, 
while s stands for the species type. Clearly, è = 1,...,9 and s = 1,..,k. As antici- 


pated, we shall be concerned with migration between adjacent cells a process which is 
encapsulated in the following relations: 


Nip & Xx HI 
Ei+xXÎ 2% Bi + Xi 


where j and 7’ represent nearby sites and E stands for the empty case. Finally, we shall 
be accommodating for the diffusion (in/out) through the region delimited by the shell. 
Such a diffusion occur either towards the inside of the cell or the outside. These effects 
are specified by: 


Bs,out 


i Bs,in i 
E° — Xi 


The allowed reactions are also depicted in Fig. 5.6 Introducing ni to label the popula- 
tion amount of species s living in cell 7, one can characterize the status of the system as 
a vector n = (n°, n°,...,n°) where n° = (ni, ni, ..., ni). With analogy of the previous 
section, we simplify the model assuming that all the rates rs, 4, 0s.in and Bs out are the 
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Figure 5.6: We consider a tiny, inner, shell adjacent to the external membrane. Such a shell is par- 
titioned in Q micro-cells (dashed yellow sectors). In every cell, autocatalytic reactions do occur of 
the Togashi and Kaneko type. As an additional ingredient, migration between nearby cells is here al- 
lowed. An effective diffusion mechanism is also invoked to regulate the mass exchange between every 
micro-cell and the interior bulk (resp. outside environment). 
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Figure 5.7: Time evolution of species 1 (red line) and species 2 (blue line) in cell number 1. Parameters 
used for the stochastic simulation are a = 0.1, r = 10, Bin = 10/64, Bout = 10/64, Q = 16, k = 4, 
N = 1000. 


same for all the chemical species, and we will drop the subscript s on the constants.The 
rates of transition from one state to the other as controlled by the above chemical equa- 
tions, are listed below (adopting the standard convention). 

Transitions stemming from the autocatalytic cycles: 


Tn, - 1, 541 + lng ns) = On ae 
For the outward/inward incoming/escape: 
i i Bout ni 
T(ni-1jni) = È 
(mi-in) = Pet 
i è Din Sa 
T(ni+1jni) = oy) 


where use has been made of the equation of conservation of mass (all species amount, 
including the empties, should sum up to N). Finally, the transition rates associated to 
the cell-to-cell migration are: 


- sr ; I a n? k n? 

T(n? -1 n? +1\n?2,n2 = S om 
(nz » Ms + [nl n? ) zQ mal 2 =) 

j i ee a ni k ni 

|P J 1 J. =] J J = SS ~ L= n 
(nz + DELA ai ) zQ N ( 2 N ) 


where z stands for the number of nearest neighbours. 

Also in this case, we can use the transition probabilities to simulate the model with 
the Gillespie algorithm. Fig. 5.7 shows a typical result of the stochastic simulation rel- 
ative to cell number 1 for a system populated by four species. Red line represents the 
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time evolution of specie 1, while the blue line refers to the second specie. With analogy 
with the non-spatial model, we observe that even-index species are asynchronized with 
those with odd-index. 

The evolution of the system is governed by the master equation which can be cast in 
the form: 


© P(n,t)= A+B+0 


where the A, B and C refer to the processes isolated above and explicitly read: 
Q k 
A= Y > [Ti ni ni + 1,05, -1) P(r) +1,n3,,-1,4) 
j=1 s=1 


-T(n} = bia + sana t)| (5.37) 


Q k 
B = LL D [Toin ni +n -DPn 10) 
j=l s=1 j'ej 


-T (ni -1,nî + 1[ni,ni )P(ni,nÎî,t) 
+T(ni,nî|ni-1,nÎî + 1)P(ni-1,nî +1,t) 
-T(ni + ln? - ini, ni )P(ni, ni D] (5.38) 


Qk 
C = YY [Tint +1) P(ni + 1,8) - T(ni - 1\nd) P(r, t) 
j=l s=1 


+T(nj|n} - 1) P(ni - 1,6) - T(ni + Ind) P(ni, t)| (5.39) 
Introducing the shift operators as 


et f(....mt,...)= fC. ni 1...) (5.40) 


the master equation takes the form 


d 2o i 
a) = 2 lG: 8 E #13; -1)T(n} - 1 nÍ +I, ni) 
SP | 
Q k P sagt 
+YL Ley DT - Lf! + ni, nf) 
jl s=1 j'ej 
xP(ni, nit) + 
+(€5 E3 7 — DT (ni + 1,nî -1n n? )P(ni, n! w] 


4 2 D [E3 DT - Unt P(24,2) 


+(Ezy -IT (nå + 1nd) P(n, 6)| (5.41) 
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We now proceed according to the van Kampen prescriptions, splitting every discrete 
variable n? into two components 


nj = Noi +N? 


where œf is the deterministic part defined as (¢”)/N, while €! is the stochastic contri- 
bution. Within this change of variables, the step operators (5.40) have a simple form for 
large N, namely 
2 
ey Sal Ni anila 
' Ogs Og)? 
= 1+ N30 + NO +... 


In this way, defining the new probability distribution II(€2,¢) = P(ni, t), the left hand 
side of the master equation (5.41) may be written as 


L Po, t)= Sige N?2 N°), DEE a 


The right hand side of the master equation can be also expressed as function of the 
expanded operators and of the stochastic variables, II, ©} and 1. We start by noticing 
that 


nio 


(eijn) 2 N (3g -3g ) +5] Hog - de, )] 


(ESES -1) 2 


1? 
NI 
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D 
| 
n° 
n 
+ 
l 
z 
l 
NI! 
— 
D 
n 
3 
LL 


12 


(€5,3€5,9 — 1) (04 =o i) F 3[N-*(0, =0 :)] 


where the operators Lı, and Lo, respectively read 
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In = (04-9%.,,) 


and 


DI 
DS 
Il 
~ 
D 
| 
D 
Z 


a) 
X 
| 
~ 
D 
D 
Z 


90 Francesca Di Patti 


In addition 
(e,,-1) = -N ða + N82 
sj 2 309 E 
839 


1 1 
+ DE: Ti . = Ape a2. 
(ef ,-1) >= Og + N Oe 


For each of the three terms (5.37)-(5.39), namely A, B, C, we can identify the associated 
N-!/? and N°! contributions. 


5.2.2 Right hand side of the master equation: N~'/? terms 


The leading term relative to A mentioned above reads: 
Ayan = 5 DRC 
Recalling the definition of La one obtains 
Ay-1/2 = “LLG CLARA dei JI 


but 


DL #9, TES En 


j s=1 j s=2 


and since k + 1 > 1 one finally gets 
Ay-p = > La (Pigi, -pga ill 
We now turn to considering the contribution of term B. 


By = -5 ALL Li(#0- Lol) +o- > on) a (5.42) 


Ss Jj ‘ej m 


which, recalling the definition of LZ, j» implies 


ee HUEY aaa- Don) -dA-Ldh, ))n (5.43) 


s j'ej 


where the factor 2 comes from exchanging j’ > j in the term associated in (5.42) to 
o) cs We then proceed by adding on the right hand side two terms which sum to zero, 


namely 


0= $3 >) bm PD Pm 
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a formal step which allows us to rewrite equation (5.43) in terms of differences between 
homologous quantities evaluated at the adjacent site 


Bien = 2° EE La[t- +0 - DEA 


s j'ej 


+i Y(di, - o] II 


m 


Introducing the discrete Laplacian 


Af; = : Xr- fi), (5.44) 


ej 


eventually yields 
Bn =-25)] Aoi. Foi.) +64 Ao}, |0g1 
Ji. 38 m 
Let us now consider the contribution relative to C for which one immediately gets 


Cram = DE (fato: fea- yolan 


Retaining order N~'/? terms in the development of master equation we obtain the 


mean-field equation for species s in cell j 


ae = pad i Eeh) tAE Adi) 


m 
-Bg 


e 2 phh) (5.45) 


Here time is re-scaled as 7 = +. Notice that by setting a to zero we immediately recover 
the mean-field equations (5.11) for the non-spatial model. 

To find the homogeneous equilibrium solution, we have to impose that no gradient 
in the concentration is allowed between nearby micro-cells. Hence the Laplacian con- 
tribution in equation (5.45) can be set to zero. The equilibrium point of the dynamics is 
therefore 

Bin — 


v= kBin + Bout F Bout 
for any s = 1,..., k and j =1,...,Q. 
5.2.3. Right hand side of the master equation: N! terms 


It should be remarked that for any of the three contributions entering the master 
equation, namely A, B and C, two types of terms are to be considered 
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e The terms werefi rst derivatives enter. These are expressed as a function of the 


operator L,,, so function of thefi rst derivatives in the fluctuation variables. These 
contributions are here called Al,_;, Bi: and Ci, 


e The terms where second derivatives enter. These are expressed as a function of the 


operator Logi so function of the second derivatives in the fluctuation variables. 
These contribution are termed A*,_,, BY, and Cy; 


Let us start by evaluating the above contributions. 


Evaluating A},_,, BY, and Ci, 


Consider the autocatalytic reaction. We have 
Aya=DLq 5 [Eis o4ele + #2162) JO 
and, recalling the definition of Lig we obtain 
Arab Dg elhia + Hat) — Os, (Hn + oF 8) 
Playing with the fe m like we did before eventually results in 


Abe = DEG (0 6 (Gas CAE 


The term B1, corresponds to 


Bra = SELD ilal- e(z) 


s j'ej 


ne (26)c (200 


which after making it explicit La reads 


By- = [DEE Oe EA -EL bn + GE 


vol D+F Loh) W-LY E g a BL 
“EY bin + & ef +g) ye. Eroh) rl (5.46) 


Changing the index j’ to j in the second combination of nested sum in equation (5.46) 
we can re-write the expression for Bi,_, as: 


BL, - = pPPLACE ef -GLE-ELH+E LE, 


s j'ej 


+ Y dh) nl 
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We can then re-write the above equation in a completely equivalent form by adding in 
each of the sum a pair of fictitious quantities, each summing up to zero 


Ss 2 On 
di > ĉi,- pi DE E: 


m 


0 


0 


Re-ordering the terms so to bring in evidence differences of species (both mean-field 
and fluctuations terms) between consecutive micro-cells and recalling the definition 
(5.44) of discrete Laplacian yield 


bue (Ae +6 Ads, - Ads Deh, + DAG, 
-ag oi,) 


Finally we focus on C},_; and immediately get 


Che =D Ed (Bel + ae ye) 


m 


Evaluating A?,_,, BX, and C%_, 


The contribution due to the second derivative in A is 
AȘ- LR Zi FLL 5 LT, (6 i IL 


which equivalently reads 


a. E ; 0? 0° _ 0° 
sy 24 Luz Ia tee | 


Then as concerns By, 


Q 
Bia = ali. 


sarigi ze) za) 


m 


VA 
(as “or e)? 
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It should be noted that the following relation a 


0? 
YE [a sal = ary UE)? Og" -agoe Jare] (5.47) 


where J<jj'> is equal to 1 if j’ and j are nearest neighbors, zero otherwise. This is a 
technicality which enables us to cast in a compact form the coefficients of the matrix 
B defined below. It should be noted that the factor 4 in the coefficients of matrix B 
originates from the factor 2 in equation (5.47). 

Finally the C3._, term is 


RARR 3 arl det) 


ente) ila 


5.2.4 The Fokker Planck equation 


Equating the orders N! in the Master equation leads to the following Fokker Planck 
equation which governs the probability distribution function of fluctuations: 


on FLL a +S 2 TE por 


m 


de dei vale tal 
The matrices A and B are defined on the basis of the above calculations. In particular we 
shall concentrate on the analysis of fluctuations around the fixed point. This is obtained 
by imposing %1 = ¢* for each choice of s and j. In the following, the structures of the 
matrices A and B are discussed. Their elements are expressed as an explicit function of 
the chemical parameters specifying the model at hand. 

The matrix A has dimension kQ x kQ and it can be explicitly written as 


Ao Ci 0 0 Ci 

Ci Ao Ci 0 0 

0 Cy Ao Cy 0 

A= 0 0 Cy Ao 0 
0 0 0 O Saxe Ci 

Ci 0 0 0 Ao 

where Ao is a circulant k x k matrix 

Ao ay, ag ag ... a3 

a3 ag Aa, ag ... Q2 

a2 a3 do ay a2 

Ao = a2 a2 a3 do «ee a2 
a2 Q2 Q2 Q2 ai 
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and C; is also k x k and takes the form 


Co Ci 
Ci Co 
Cr Ci 
Ci Ci 
Cy CI 
Ci Ci 


Q 


Co 
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The elements of such matrices follow from the calculation developed above and read 


— 2a Bin Bout 
2a, Ta Bin 
Oy A 
a a -24g Bin 
o Q 
2A Ta Bin 
OS al o 
2 
o = Gl-(k-1)01 
am $ 
Cy = 79° 
The coefficient of matrix B are 
259° dd; if |q - s| = 
[220*¢* + p (1- A 
Bij qs = eda (1 ko") + Sed] diy ifq=s 
-42 (1 - ko”) Jejjt> 
0 otherwise 


More specifically matrix B is kQ x kQ and has the following structure: 


Bo 
By 


Bı 0 
Bo Bi 
Bı Bo 
0 Bb, 
0 0 
0 0 
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where Bo is a k x k matrix defined as: 


bo bi 0 0 0 bi 
bi bo b 0 0 0 
0 b bo b 0 0 
Bo = 0 0 b bo 0 0 
0 0 0 0 be! By 
b 0 0 0 bı bo 
and B, is also a k x k matrix: 
b 0 0 0 
0 b 0 0 
Bı = 0 0 b 0 
0 0 0 bo 
The coefficients by, b1, b2 read: 
do = 2h prot Bpa — apt) + Pini ng) tg 
bi = 2300 
a 
bo = -4—o*(1-kd* 


Having defined the matrices A and B, we are in the position to obtain an analytical 
expression for the power spectrum for the fluctuations. This is achieved via a straight- 
forward procedure which is essentially analogous to that reported in section 5.1.1. Just 
a few technical points need to be carefully handled, as clarified in [84]. Here, we are 
solely concerned with presenting a preliminary gallery of results which testifies on the 
predictive ability of our analysis. 

In Fig. 5.8 the (two dimensional) analytical and numerical power spectra are dis- 
played, for the case where k = 4. Clearly, the power spectra are now function of the 
time-frequency w and space-frequency x. The migration parameter a is set to zero so 
that the system is practically composed by Q, independent replica of the Togashi and 
Kaneko setting. The other parameters are set to the same values as adopted in Fig. 5.3. 
As expected, a peak at k = 4 is observed, and no spatial modulation recorded. Direct 
simulations agree with the theoretically predicted profile. We then turn on the migra- 
tion effect and set a = 0.1, see Fig. 5.9a. The system still exhibits the peak for w = 4, but 
now a decay in k is predicted to occur. This feature is then observed in the simulations 
(see Fig. 5.9b), again pointing to the validity of the perturbative development, and sug- 
gesting that the system organizes collective modes on the large wavelength scale. We 
notice that a similar effect is also reported by Lugo and McKane in [86] for their spatial 
predator-prey model. 

Preliminary simulations (not reported here) where the species are differentiated 
with respect to their chemical activity (different rates) shows a large zoology of allowed 
spatio-temporal patterns, including the existence of isolated peaks in the power spectra, 
see [84] for an extended discussion. 
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(a) (b) 


Figure 5.8: Numerically (panel (a)) and analytically (panel (b)) obtained power spectrum for a system 
with k = 4. Parameters used for the simulations are r = 10, a = 0, Bin = 10/64, Bout = 10/64, Q = 16, 
N = 5000.Th e numerical power spectrum is obtained from 100 runs of the stochastic simulation. 


As a final remark, we wish to stress that the emergence of regular spatio-temporal 
oscillations in the concentration amount might have important implication for the dy- 
namics of the so-called protocells. These are small-cell like, living, units, which can 
self-assemble, develop and replicate. Autocatalytic chemical reactions might have oc- 
curred inside those primordial containers back to the origin of life, and could have hy- 
pothetically contributed to mediate the minimal mechanisms involved in duplication. 
These aspects are little understood, despite the fact that lipidic vesicles (i.e. the candidate 
representative of the protocells’ family) are nowadays heavily studied in laboratories. In 
a future work we intend to explore the possibility that the division process is initiated 
by a Turing instability, i.e. following from a purely spatial effect of the type investigated 
in this chapter. 


Figure 5.9: Numerically (panel (a)) and analytically (panel (b)) obtained power spectrum for a system 
with k = 4. Parameters used for the simulations are r = 10, a = 0.1, Gin = 10/64, Bout = 10/64, 
Q = 16, N = 5000.Th e numerical power spectrum is obtained from 100 runs of the stochastic 
simulation. 


Conclusions 


Microscopic systems, such as the intracellular environment, are characterized by 
pronounced biochemical oscillations that the standard deterministic approach fails to 
capture. This observation has stimulated the research of novel theoretical frameworks 
to interpret the large class of intriguing scenarios characterizing these systems. In par- 
ticular, the discrete nature of the system, which can be microscopically segmented in 
independent entities, can yield to an effective stochastic resonance, resulting in ampli- 
fied regular cycles. Finite size contributions do matter and may hence cause complex 
dynamical patterns for the populations under inspection. 


In this thesis we have investigated this crucial aspect with reference to a selected 
gallery of problems. First, focusing on pain perception, we have developed a discrete 
dynamical framework to study the molecular processes which are activated in response 
to an external harming stimulus. These cascade of reactions ultimately triggers the pain 
sensation. The problem has been analyzed, paying special attention to the associated 
finite-size effects. Through the van Kampen perturbative theory, we are able to recover 
an exact analytic description of the fluctuations. We were in particular interested in elu- 
cidating the crucial interplay between the administered drug molecules, which express 
their analgesic function chasing the target receptors, and other chemical elements freely 
diffusing in the bloodstream. The latter can substantially reduce the anaesthetic effect, 
by hindering the available binding sites. Similarly, drug molecules can be turned into 
inactive species following binary encounters. The mechanisms here postulated were 
formally coded via chemical reactions and define a consistent stochastic scheme. Nu- 
merical simulations displayed macroscopic oscillations in the concentration amount: 
the number of bound receptors changed cyclically in time, a trend which we assumed 
to induce an analogous modulation for the experienced perception of pain. It is impor- 
tant to remark that the amplification process here discussed stems from the underlying 
stochasticity, which is resonant with the natural frequencies of the system. Oscillations 
arise hence spontaneously, driven by the noise which is intrinsic to the system and with- 
out invoking any ad hoc couplings among the molecular agents participating to the 
dynamics. Our findings suggest the existence of a simple, though general, molecular 
mechanism responsible for the emergence of cyclic behaviors in response to analgesic 
treatments. 


An extension of the model included drug-metabolite interactions where different 
populations compete for the same target receptors, so to induce analgesia.Th e system 
has been analyzed, focusing first on the mean-field dynamics (N + 00) which is gov- 
erned by a set of coupled ordinary differential equations for the species amount. The 
fixed points have been studied together with their associated stability properties. The 
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chemical parameters are found to control the asymptotic regime determining the ef- 
fectiveness of the administered therapy. More interestingly, in the transient dynamics a 
lethargic phase is registered, where the number of bound receptors appear to have stabi- 
lized to a quota sensitive to initial condition, before reaching their equilibrium solution. 
Moreover, fluctuations have been also analyzed via the van Kampen technique. It is here 
speculated that, in particular cases, they might influence the degree of experienced pain, 
which could hence change over time. 


Another context where intrinsic noise may have important effects on the dynamics 
is that of the auto-catalytic networks. These kind of reactions are central in many differ- 
ent contexts and play an important role in intracellular biochemical reaction schemes. 
In this latter scenario, species are confined in a closed volume, delimited by the cellu- 
lar membrane. Low concentration can occasionally develop resulting from the com- 
plex mutual interaction between microscopic actors. Under such conditions, fluctua- 
tions matter and the effects of the intrinsic discreteness need to be properly accounted 
for. In other words, continuous kinetic equations prove inadequate, finite size correc- 
tions becoming significant. These aspects were numerically substantiated by Togashi 
and Kaneko [77, 78] within the framework of a simplified system of k coupled autocat- 
alytic reactions. In this thesis we elaborate on this concept by studying analytically the 
associated master equation via a systematic expansion in power of N~'/?, where N is 
the system size. To leading order, the mean-field rate equations are recovered, while 
higher order corrections enable us to explain the large amplitude of the oscillations as 
detected in direct simulations. Importantly, the calculation applies to arbitrary values of 
k. For k = 4a peak in the power spectrum is found, while for k = 8 two peaks develop. 
To the best of our knowledge, this is thefi rst time that a double-peaked power spectrum 
has been predicted to emerge as a resonant effect, within a van Kampen type of anal- 
ysis. In both cases, theory and simulations agree well thus confirming the importance 
of finite N contributions. Later on we have extended our work by taking spatial effects 
into account. Also in this case, we detected spatio-temporal oscillations in the species 
concentration, which are again driven by the discreteness of the system components. 


In addition tofi nite-size effects problems, a short part of this thesis has been devoted 
to describe new method for analyzing experimental data relative to drugs kinetics and 
to microarray experiments. 


First, we have proposed a mathematical model for the kinetics of tramadol. This 
novel theoretical framework could result in an objective criterion on how to adjust the 
assigned dose, depending on the genetic polymorphisms of CYP2D6.Th e model de- 
scribes the coupled dynamics of tramadol and the metabolite O-desmethyl-tramadol. 
The effect of diffusion of the drug in the blood is here accounted for and we further 
hypothesize the existence of a time delay in the process of chemical translation from 
tramadol into metabolites. The system of coupled differential equations is solved nu- 
merically and the free parameters adjusted so to interpolate the experimental time se- 
ries for the intravenous injection setting.Th eoretical curves are shown to reproduce 
correctly the experimental profiles obtained from clinical trials. This enables in turn 
to extract an estimate of the metabolization rate. A difference in metabolization rate 
between CYP2D6 poor and extensive metabolizers is also found, and the stereoselec- 
tivity in the O-demethylation of tramadol highlighted. Our results allow one to quantify 
the dose of (+)-tramadol (resp. (—)-tramadol) administered to poor or extensive me- 
tabolizers, if the same effect is sought. The latter is here quantified through the blood 
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concentration of (+)-metabolites (resp. (-)-metabolites). 

As for the microarray experiments, we have developed a method for measuring the 
distance among records based on the correlations among available data, with the aim 
of estimating the missing values. Finally, we have reviewed the Van Dongen algorithm, 
suggesting a new method for detecting different levels of communities, i.e. clusters of 
homogeneous objects with respect to a predefined norm. 


Chapter A 
A technique to simulate the full stochastic process: The Gillespie 
algorithm 


In 1976 Daniel T. Gillespie proposed an algorithm to exactly simulate the stochastic 
dynamics of chemical reactions [79, 82]. To describe this method, let us consider a 
volume V which contains molecules of N chemically active species S; for i = 1,..., N, 
and denote by X; the current number of molecules of species S; in V. The molecules 
interact according to M chemical reactions R, for u = 1,...,M, each characterized 
by a reaction parameter c,,. The quantity c, dt represents the first order approximation’ 
of the average probability that a particular combination of R, reactant molecules will 
react accordingly in the next time interval dt, as it follows by a chemical kinetics theory 
derived into details in the original paper [79]. 

To illustrate the relationship between c,, and the more familiar “reaction rate con- 
stant” k,, used in the deterministic formulation of chemical kinetics, let us consider 
the reaction Sı + S2 — 253. In this case, X1 Xə - c dt is the probability that the 
reaction will occur inside V in the next time interval dt, where XX represents the 
distinct combinations of reactant molecules in V. Averaging over a set of stochastically 
identical system, and diving by V, we obtain the average reaction rate per unit time 
(X,X2)c,,/V or, in terms of molecular concentrations x; = X;/V, (x1 %2)c,,V. If we 
divide this latter quantity by the product of the average densities of the reactants, we 
obtain the expression for k,,, namely 


(t102)CyV 
(71)(22) 


In the deterministic formulation the average of a product is equivalent to the product 
of the averages, thus (1122) =( x1)(x2) and (A.1) simplifies to 


k, = (A.1) 


ky = Vey 


The factor V in this relation is due to the type of reaction considered. In reactions with 
only one reactant molecule, indeed, the factor V would be absent, while in those with 
three reactants, a V° would instead appear. 

The aim of the method is to simulate the time evolution of the N variables X; know- 
ing their initial values X;(0), the M reactions R,, and the associated reaction parame- 


! More precisely, the first order in ôt means that the average probability is c,ét + o(ôt) with 
lims:—0 0(6t)/dt = 0. 
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Reaction h, 

Sj — reaction products Xj 

Sj + Sk — reaction products X; Xx 

25; — reaction products X;(X;-1)/2 


Si + Sj + Sk — reaction products X,X;Xx 
Sj + 2S), — reaction products X;Xx(X;,-1)/2 
3S; — reaction products X;(X;-1)(X; - 2)/6 


Table A.1: A selection of possible reactions and their corresponding state variables h,,. 


ters c,,. The standard stochastic approach to this problem focuses on the master equa- 
tion, namely the time evolution of the probability function P(X,,..., Xy;t) to have 
X; molecules of S; (for è = 1,..., N) at time t. In most cases this approach turns out to 
be intractable, both analytically and numerically. To overcome this problem, Gillespie 
proposed a method based on what he called the reaction probability density function 
P(7, 1). He defined this quantity as the probability at time # that the next reaction in 
V will occur in the time interval (t + 7,t + 7 + 67) and that the selected reaction was 
of the type Rp. 


Thefi rst step for deriving an analytical expression for P(r, u), consists in associ- 
ating to every chemical reaction, a state variable h, defined as the distinct molecular 
reactant combinations for reaction R, within the volume V at time t. Table A shows 
the state variables for a selection of reactions. In this way h,,c,,dt is the probability, to 
first order in ôt, that an R, reaction occurs in V, in the next time interval ôt. 


The second step requires decomposing P(r, u) as 
P(t, w)dt = Po(t)- hycydT (A.2) 


where P(T) is the probability at time t that no reaction will occur in the time interval 
(t,t + 7), and h,,c,,d7 is the probability that an R, reaction will occur in the next 
differential time interval (t + 7,t+7+ dr). 

To calculate P)(7) one can divide the interval (t,t + 7) in K subintervals of equal 
length e = T/K. In each subinterval, the probability that none of the reactions occurs 
is given by 
1 M 
[1-A,cye+o(e)]=1- E hycye+o(e) (A.3) 


1 v=1 


an 


Vv 


In this way, Po(7) is just the product of K times equation (A.3) 


P(T) | T= > hucre + 7) 


v=1 


M = ts 
= |1- hve +0(K*) 
piane 
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This relation holds for any K > 1, and therefore it is true for infinitely large value of K: 


K 


K-00 


oh M -1 È 
Boa da È Li hyeyT + 0(K = 


M 
= exp |- X hesr (A.4) 
v=1 


Putting together equation (A.2) with equation (A.4) one obtains the exact expression 
for the probability density function 


M 
P(T, u) = hC, exp |- di hes (A.5) 
v=1 


for 0 < T < œ and 1 < u < M with 7 € R and u € N. 

Before moving to the description of the algorithm, we recall the main ideas of a 
Monte Carlo method. This latter constitutes a crucial step in the Gillespie implemen- 
tation, providing a method to generate two random numbers 7 (real) and u (integer) 
according to the joint probability density function in (A.5). The trick consists in splitting 
the probability density function P(r, 4) into the product of two one-variable proba- 
bility density functions.Th is procedure is called conditioning and leads to 


P(r, u) = Pi(7)- P2(ul7) (A.6) 


where P(7)dr is the probability that the next reaction will occur between times t + 7 
and t+7 + dr, and P;(p|r) is the probability that the next reaction will be an R, type, 
given that it happens at time t + 7. Invoking the addition theorem for probabilities, 
P,(T)dr is obtained by summing P(T, 4) d7 over all ju, and thus 


M 
Pi(r) = XY P(T, u) 

pal 
Putting this into (A.6) and solving for P2(,u|T) it gives 


P2(ulr) = P(r, p)/ DE) 


Substituting P(7, p) with (A.5) in the previous two equations, we obtain 


_f aexp[-at] for0<T <0 
Pur) = { 0 otherwise (A.7) 
and z 
|) aul- fov=1,..., M 
P,(ulr) = | 0 otherwise (39) 
where È 
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with 
ap =hpernw foru=1,..., M (A.9) 


In this way the problem of finding two random numbers according to P(T, 4) may 
be recast as the problem of drawing a real random number from the P, distribution, 
and an integer random number according to P}. 

Let us first focus on the former case. We wish to generate a real number x accord- 
ing to a probability density function P(x). The corresponding probability distribution 
function 


F(a») = / si P(a)dx (A.10) 


quantifies the probability that x will be less than xo. The inversion method for gen- 
erating a random value x according to P(x) is to draw a random number r from the 
uniform distribution in the unit interval, and than take 


x=F3(r) (A.11) 


To prove that this procedure is correct, we have to show that the probability that the x 
value so generated will lie between x’ and x’ + dx’, is P(x')dx'. By construction, this is 
equivalent to calculating the probability that r will lie between F'(x’) and F(a’ + dz’). 
Since r is a random number drawn from the uniform distribution in the unit interval, 
this probability is just the length of the interval | F(x’), F(a’ + dx')], namely F(x’ + 
dx’) — F(x’) = F'(x')da'. Applying the definition (A.10), we get 


F(x' + da')- F(x') = F'(x')dx' = P(x')da' 


and this prove that the probability density function for the random number x generated 
according to (A.11) is indeed P(x). 

For the specific case at hand, we wish to generate a random number 7 according 
to the probability density function (A.7). In this case F(7) = 1- exp[-ar]. Putting 
F(T) = r and inverting the function F, we obtain 


T= ~In(=) (A.12) 
a r 
where, for simplicity, we have replaced the random variable 1 — r by the statistically 
equivalent random variable r. 

We have seen how to generate a random number according to a specific probability 
density distribution for a continuous variable. Now we consider the discrete case and 
we look for a method which enables us to obtain a random integer è according to the 
probability density function P(j), where now P(j) is the probability that è = j. The 
corresponding distribution function F (i) is defined by 


FO- $ PO) 


and F(io) represents the probability that 7 < iọ. With analogy to the continuous case, 
the inversion method consists in drawing a random number r from the uniform distri- 
bution in the unit interval and take for è that value which satisfies 


F(i-1) <r< F(i) (A.13) 
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To show that the procedure is correct also in this case, we use the fact the resulting 
integer i will equal j is equivalent to the probability that r will lie between F(j-1) and 
F(j). So we have 


F()-FG-1)= Xx P(k) - > P(k) = P(j) 


This proves that P(7) is indeed the probability density function for the random integer 
i generated according to (A.13). 

As an example, we consider again our specific case, and make it explicit the expres- 
sion of the random integer u with respect to the density function (A.8). Applying (A.13) 
we see that we have to select the integer u so that 


pel H 
X Po(v|r) << Y P(r) 
v=1 v=1 


or 
pel M H 
> Gxt > as’ a (A.14) 
v=1 v=1 v=1 


Now we have all the ingredients to describe the details of the simulation methods. 
The steps of the algorithm are the following: 


Step0 Assign values to the M reaction constants c1, . . . , cy and initialize the N molec- 
ular population numbers X1,..., Xy. Set the time variable t = 0, and specify a 
stopping time tstop- 


Stepl Calculate the quantities a, = h,c, for v = 1,..., M for the current molecular 
population numbers, and the quantity ao = YM, a,. 


Step2 Use the Monte Carlo technique to generate a random pair (7, 4.) according to 
(A.12) and (A.14). 


Step3 According to the numbers 7 and u generated in the previous step, advance time 
by 7 (t = t +T) and update the values of X; for every species involved in reaction 


i 
Step4 Ift < tstop go to stepl , otherwise terminate the calculation. 


It is important to stress that the time series generated with this algorithm recover the 
exact probability distribution function given by the master equation. It can be shown, 
in fact, that the two approaches, the master equation and the Gillespie's method, are 
equivalent at thefi rst order approximation grounded on the kinetic theory argument 
mentioned before. 
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